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EXECUTIVE  SUMMARY 


The  research  conducted  under  this  grant  focused  on  the  design  and  analysis  of  new  algorithms  within  the 
generalized  hill  climbing  algorithm  framework  for  address  intractable  discrete  optimization  problems. 
The  major  technical  accomplishments  achieved  include 

i)  the  introduction  of  the  simultaneous  generalized  hill  climbing  (SGHC)  algorithm  framework  for 
addressing  sets  of  related  intractable  discrete  optimization  problems, 

a)  the  formulation  of  sufficient  conditions  under  which  SGHC  algorithms  converge, 

Hi)  the  formulation  of  finite-time  performance  measure  for  generalized  hill  climbing  (GHC) 
algorithms. 

iv)  the  application  of  such  measures  to  a  static  simulated  annealing  algorithm,  which  demonstrated 
how  the  performance  of  such  algorithms  can  be  predicted. 

v)  the  application  of  GHC  algorithms  to  a  construction  site-leveling  problem. 

Several  other  problems  were  studied  during  the  grant  period,  including  estimation  procedures  for 
discrete  event  simulation  problems,  aviation  security  system  design  and  optimization  issues,  and  a 
pediatric  vaccine  formulary  design  problem.  All  the  accomplishments  are  documented  in  several 
archival  journal  articles  and  conference  proceedings.  In  addition,  many  of  the  results  have  been 
presented  at  national  and  international  conferences,  and  have  won  awards  for  their  contribution. 

Three  Ph.D.  dissertations  were  completed  during  the  period  of  the  grant.  Colonel  (Sel.)  Darrall 
Henderson  successfully  defended  and  submitted  his  Ph.D.  dissertation  "Assessing  the  Finite-Time 
Performance  of  Local  Search  Algorithms"  in  April  2001.  Dr.  Tevfik  Aytemiz  successfully  defended  and 
submitted  his  Ph.D.  dissertation  "  A  Probabilistic  Study  of  3-SATISFIABILITY"  in  July  2001.  Dr. 
Derek  Armstrong  successfully  defended  and  submitted  his  Ph.D.  dissertation  "A  Local  Search  Algorithm 
Approach  to  Analyzing  the  Complexity  of  Discrete  Optimization  Problems"  in  May  2002. 
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1.  Generalized  Hill  Climbing  Algorithms 

Generalized  hill  climbing  algorithms  (Jacobson  et  al.  1998)  provide  a  well-defined  framework  to  model 
algorithms  for  intractable  discrete  optimization  problems.  Generalized  hill  climbing  algorithms  allow 
inferior  solutions  to  be  visited  enroute  to  optimal  solutions.  This  hill  climbing  capability  is  the  basis  for 
the  search  strategy's  name.  To  formally  describe  the  generalized  hill  climbing  algorithm  framework, 
several  definitions  are  needed. 

For  a  discrete  optimization  problem,  define  the  solution  space,  Q,  as  a  finite  set  of  all  possible 
solutions.  Define  an  objective  function  f:  Q  ->  [0,+oo]  that  assigns  a  non-negative  value  to  each  element 
of  the  solution  space.  An  important  component  of  GHC  algorithms  is  the  neighborhood  function,  t]:  Q 
2^,  where  ri(oo)  c  for  all  co  e  Q.  Solutions  in  a  neighborhood  are  generated  uniformly  at  each 
iteration  of  a  GHC  algorithm  execution  if,  for  all  ooeQ,  with  co'  €  r|((o), 

P{(o'  is  selected  as  the  neighbor  of  co  at  a  given  iteration  of  a  GHC  algorithm}  =  1  /  |  ti((o)  | . 

Unless  otherwise  stated,  assume  that  the  neighbors  are  generated  uniformly  at  each  iteration  of  a  GHC 
algorithm.  The  GHC  algorithm  is  described  in  pseudo-code  form  (Jacobson  et  al.  1998): 

Set  the  outer  loop  counter  bound  K  and  the  inner  loop  counter  bounds  N(k),  k  =  1,2,. .  .,K 
Define  a  set  of  hill  climbing  (random)  variables  Rki  Q  x  Q  u  {-oo,-i-oo},  k=  1,2,...,K 
Set  the  iteration  indices  i  =  0,  k  =  n  =  1 
Select  an  initial  solution  co(0)  e  Q. 

Repeat  while  k  <  K 
Repeat  while  n  <  N(k) 

Generate  a  solution  co  e  ri(co(i))  and  calculate  5(co(i),co)  =  f(co)  -  f(co(i)) 

If  Rk(oo(i),co)  >  5(co(i),ci)),  then  co(i-l-l)  <—  co  (accept  improving  or  hill  climbing  moves) 

If  Rk(co(i),co)  <  5(co(i),co),  then  co(i+l)  <-  co(i)  (reject  hill  climbing  moves) 

n  •<—  n+1,  i  <-  i+1 
C/«ri/n  =  N(k) 
n<—  l,k<— k+1 
Until  k  =  K 

Note  that  the  outer  and  inner  loop  bounds,  K  and  N(k),  k  =  1,2,. .  .,K,  respectively,  may  all  be  fixed,  or  K 
can  be  fixed  and  the  N(k),  k  =  1,2,..  .,K,  are  random  variables  with  values  determined  by  the  solution  at 
the  end  of  each  set  of  iimer  loop  iterations  satisfying  some  property  (e.g.,  the  solution  is  a  local  optima). 
Moreover,  assume  that  the  hill  climbing  random  variables  have  finite  means  and  finite  variances  for  all  k 
and  for  all  possible  pairs  of  elements  in  the  solution  space  (i.e.,  E[Rk(co(i),co)]  <  +oo  and  Var[Rk(co(i),co)] 

K 

<  +CO  for  all  co(i)  e  Q,  coer|(co(i)),  and  for  all  k  =  1,2,. .  .,K,  i  =  1,2,. .  .,I  =  X  N(k)). 

i=l 

The  neighborhood  function  establishes  relationships  between  the  solutions  in  the  solution  space, 
hence  allows  the  solution  space  to  be  traversed  or  searched  by  moving  between  solutions.  To  ensure  that 
the  solution  space  is  not  fragmented,  assume  that  all  the  solutions  in  the  solution  space  (with 
neighborhood  function  t|)  are  reachable;  that  is,  for  all  co',co"eG,  there  exists  a  set  of  solutions  coj,  CO2,. 
cOm  e  Q  such  that  co^  g  r|(cOr.i),  r  =  l,2,...,m+l,  where  co'  =  coo  and  co"  =  cOn^-i.  Note  that  if  all  solutions  in  a 
solution  space  are  reachable,  then  the  solution  space  (with  neighborhood  function  r|)  is  said  to  be 
reachable.  The  goal  is  to  identify  a  globally  optimal  solution  co*  (i.e.,  f(co*)  <  f(co)  for  all  co  g  Q). 

2.  Simultaneous  Generalized  Hill  Climbing  Algorithms 

One  accomplishment  during  the  term  of  this  grant  has  been  the  development  of  a  mathematical 
framework  for  simultaneously  addressing  a  set  of  related  discrete  optimization  problems  using  local 
search  algorithms.  The  resulting  algorithms,  termed  simultaneous  generalized  hill  climbing  (SGHC) 
algorithms  (Vaughan  and  Jacobson  2003a),  can  be  applied  to  a  wide  variety  of  sets  of  related 
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manufacturing,  military  and  service  industry  discrete  optimization  problems.  Many  well-known  local 
search  algorithms  can  be  embedded  within  the  SGHC  algorithm  framework,  including  simulated 
annealing,  threshold  accepting,  Monte  Carlo  search,  and  pure  local  search  (among  others). 

A  SGHC  algorithm  probabilistically  moves  between  the  set  of  related  discrete  optimization  problems 
during  its  execution  according  to  a  problem  generation  probability  function.  The  problem  generation 
probability  function  is  shown  to  be  a  stochastic  process  that  satisfies  the  Markov  property.  Therefore, 
given  a  SGHC  algorithm,  movement  between  discrete  optimization  problems  can  be  modeled  as  a 
Markov  chain.  Sufficient  conditions  that  guarantee  that  this  Markov  chain  has  a  uniform  stationary 
probability  distribution  are  obtained.  Computational  results  are  presented  with  SGHC  algorithms  for  a 
set  of  traveling  salesman  problems.  For  comparison  purposes,  GHC  algorithms  are  also  applied 
individually  to  each  traveling  salesman  problem.  These  computational  results  suggest  that  optimal/near 
optimal  solutions  can  be  reached  more  effectively  using  a  SGHC  algorithm. 

It  is  common  to  encounter  several  discrete  optimization  problems  where  a  relationship  exists  between 
the  solution  spaces  of  the  individual  problems.  In  general,  these  problems  are  approached  individually. 
However,  because  of  their  similarities,  the  same  computational  tools  can  be  effectively  used  to  address 
them  simultaneously.  For  example,  the  Material  Process  Design  Branch  of  the  Air  Force  Research 
Laboratory,  Wright  Patterson  Air  Force  Base,  has  several  similar  discrete  manufacturing  process  design 
optimization  problems  under  study  (see  Jacobson  et  al.  1998,  Sullivan  and  Jacobson  2000).  These 
problems  are  difficult  to  solve  due  in  part  to  the  large  number  of  design  sequences  and  associated  input 
parameter  setting  combinations  that  exist.  Local  search  algorithms  in  the  generalized  hill  climbing  (GHC) 
algorithm  framework  were  introduced  to  address  such  manufacturing  design  problems  (Jacobson  et  al. 
1998).  Initial  results  with  GHC  algorithms  required  the  manufacturing  process  design  sequence  to  be 
fixed,  with  the  GHC  algorithm  used  to  identify  optimal  input  parameter  settings  for  each  feasible  design 
sequence  (Jacobson  et  al.  1998). 

The  SGHC  algorithm  framework  is  motivated  by  a  study  reported  in  Vaughan  et  al.  (2000)  which 
develops  a  new  neighborhood  function  that  allows  a  local  search  algorithm  (in  the  GHC  algorithm 
framework)  to  be  used  to  also  identify  the  optimal  discrete  manufacturing  process  design  sequence  among 
a  set  of  valid  design  sequences.  This  neighborhood  function  allows  the  GHC  algorithm  to  simultaneously 
optimize  over  both  the  design  sequences  and  the  input  parameters,  eliminating  the  need  to  approach  each 
design  sequence  (i.e.,  a  discrete  optimization  problem)  individually.  The  computational  results  in  Vaughan 
et  al.  (2000)  suggest  that  such  an  approach  is  feasible  and  yields  reasonable  results.  However  the 
neighborhood  function  developed  for  this  purpose  is  overly  complex  and  problem  specific. 

The  following  analysis  generalizes  the  results  reported  in  Vaughan  et  al.  (2000)  by  formally  defining  a 
class  of  sets  of  discrete  optimization  problems  where  a  relationship  exists,  similar  to  the  one  described  for 
the  manufacturing  problem.  A  set  of  discrete  optimization  problems  contained  in  this  class  is  defined  as  a 
set  oi  fundamentally  related  discrete  optimization  problems.  The  SGHC  algorithm  framework  is  used  to 
simultaneously  approach  sets  of  fundamentally  related  discrete  optimization  problems  using  GHC 
algorithms.  A  metric  between  elements  in  a  set  of  fundamentally  related  discrete  optimization  problems  is 
introduced  to  evaluate  if  and  when  it  is  advantageous  to  address  a  particular  set  of  discrete  optimization 
problems  with  a  SGHC  algorithm. 

Vaughan  et  al.  (2000)  introduces  a  new  neighborhood  function  for  simultaneously  addressing  a  set  of 
related  manufacturing  process  design  optimization  problems  using  GHC  algorithms.  This  neighborhood 
function  allows  for  simultaneous  optimization  across  the  design  sequences  and  the  controllable  input 
parameters.  The  application  of  such  optimization  algorithms  (that  simultaneously  explore  multiple 
manufacturing  process  designs)  using  computer  simulation  is  a  new  advance  in  how  optimal  manufacturing 
process  designs  can  be  efficiently  identified  (Vaughan  et  al.  2000). 

It  is  common  to  encounter  several  discrete  optimization  problems  where  a  relationship  between  the 
solution  spaces  of  the  individual  problems  exists.  For  example,  Henderson,  Vaughan  and  Jacobson  (2003) 
introduces  a  multiple  platform  search  and  rescue  optimization  problem  that  is  modeled  as  several  discrete 
optimization  problems.  They  show  that  for  this  set,  the  solution  spaces  of  the  individual  problems  overlap 
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and  are  a  set  of  fundamentally  related  discrete  optimization  problems.  Vaughan  et  al.  (2000)  describes  an 
integrated  blade  rotor  diserete  manufacturing  process  design  problem  that  illustrates  how  certain 
manufacturing  problems  can  be  modeled  as  several  discrete  optimization  problems  with  overlapping 
solution  spaces  that  satisfy  the  definition  of  fundamentally  related  discrete  optimization  problems.  Note 
that  this  paper  relaxes  the  methodology  used  to  address  the  integrated  blade  rotor  discrete  manufacturing 
process  design  problem  described  in  Vaughan  et  al.  (2000)  to  develop  a  general  mathematical  framework 
for  simultaneously  approaching  sets  of  fiindamentally  related  discrete  optimization  problems. 

To  discuss  the  class  of  sets  of  discrete  optimization  problems  for  which  SGHC  algorithms  are 
applicable,  the  following  definitions  are  needed.  Consider  a  set  of  discrete  optimization  problems  S  =  {Di, 
D2,  ...,  Dm},  where  each  discrete  optimization  problem  Dy=  (Qy,  fy)  is  fully  defined  by  a  finite  set  of 
solutions  Qy  and  a  real-valued  objective  function  fy:  Qy->i?.  A  set  of  discrete  optimization  problems  S  is 
fundamentally  related  by  a  set  Ob  =  {ci,  C2,  ...,  c„}  of  objects  if  the  solution  space  Qy  of  each  discrete 
optimization  problem  Dy  =  (Qy,  fy)  g  S  can  be  fully  defined  by  exactly  one  subset  of  Ob.  Then  for  every 
discrete  optimization  problem  Dy  =  (Qy,  fy)  e  S,  there  is  exactly  one  set  Cy  c  Ob  such  that  Cy  completely 
defines  Qy.  The  set  Cy  is  defined  to  be  ihs  fundamental  relation  set  of  Dy. 

Let  S  be  a  set  of  fundamentally  related  discrete  optimization  problems  related  by  Ob.  Consider  Dy  e  S 
where  Cy  c  Ob  is  the  fundamental  relation  set  of  Dy.  Then  Cy  can  be  represented  by  the  binary  activity 
vector  G  {0,  1}”,  c’'=  (Bi,  B2,  ...,  Bn),  where 


fl ,  if  C;  is  contained  in  Cy 
[  0,  otherwise 


SGHC  algorithms  are  developed  to  approach  sets  of  fundamentally  related  discrete  optimization 
problems.  When  two  discrete  optimization  problems,  Dy  and  Dq,  are  contained  in  a  set  of  fundamentally 
related  discrete  optimization  problems  with  respective  fundamental  relation  sets,  Cy  and  Cq,  where  [CyOCql 
/  \Ob\  is  close  to  one,  it  is  reasonable  to  conjecture  that  the  optimal/near  optimal  solutions  of  Dy  and  Dq  are 
similar.  The  following  detachment  metric  is  designed  to  determine  if  two  discrete  optimization  problems, 
in  a  set  of  fundamentally  related  discrete  optimization,  are  close  together,  hence  have  similar  solution 
spaces. 

Let  S  be  a  set  of  fundamentally  related  discrete  optimization  problems  related  by  Ob.  To  formally 
define  the  detachment  metric  p  between  discrete  optimization  problems  Dy,  Dq  g  S,  consider  the  metric 
space  <E”,  p>,  with  S  =  {0,  1},  where  the  detachment  metric  p  is  defined  on  S"  x  E”  such  that  the  distance 
between  two  discrete  optimization  problems  can  be  determined  by  considering  their  binary  activity  vectors 
,cl,  ...,cl)  G  E”and  €'•=  (Ci^,c|,  ...,c^)gZ".  ihs  detachment  metric  zs 


p(Dy,  Dq) 


c  -C^\  + 


+  ...+ 


-c‘‘ 


(Royden  1988,  p.l40).  The  detachment  metric  provides  a  way  to  measure  the  overlap  (or  lack  of  overlap) 
between  the  members  in  a  set  of  fundamentally  related  discrete  optimization  problems. 

To  apply  SGHC  algorithms,  a  neighborhood  function  with  an  associated  problem  generation 
probability  function  for  moving  between  discrete  optimization  problems  during  an  execution  of  a  SGHC 
algorithm  must  be  developed.  The  neighborhood  function  is  defined  such  that  each  discrete  optimization 
problem  has  the  entire  set  of  discrete  optimization  problems  as  neighbors.  Therefore,  whenever  this 
neighborhood  function  is  applied,  every  discrete  optimization  problem  is  a  candidate  problem  (i.e.,  has  a 
positive  probability  of  being  selected).  The  problem  generation  probability  function  determines  the 
probability  of  selecting  a  candidate  problem. 

More  formally,  define  the  neighborhood  function,  rjset:  S  — >  2®,  such  that  r|set(Dy)  =  S,  for  all  Dy  g  S. 
Define  the  problem  generation  probability  function  (k,  p(Dy,  Dq)),  such  that 

0  <  2),  P(^y>  Dq))  <  for  every  Dy  g  S,  Dq  g  Tiset(Dy), 
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where 


X  X,  (k,  p(Dy,  Dq))  =  1 ,  for  every  Dy  e  S,  Dq  g  riset(Dy) 

for  every  k  =  1,  2, . . K. 

Note  that  the  problem  generation  probability  function  (the  probability  of  selecting  a  candidate  problem, 
Dq  e  riset(Dy),  Dy  e  S)  can  be  a  function  of  both  the  outer  loop  iteration  k  =  1,2,...,K  and  the  detachment 
metric  p(Dy,  Dq). 

SGHC  algorithms  provide  a  mathematical  framework  for  addressing  several  fundamentally  related 
discrete  optimization  problems  simultaneously  using  GHC  algorithms.  SGHC  algorithms  seek  to  find 
optimal  solutions  for  sets  of  fundamentally  related  discrete  optimization  problems  by  allowing  the 
algorithm  to  probabilistically  move  between  discrete  optimization  problems.  When  a  new  discrete 
optimization  problem  is  generated,  an  initial  solution  for  this  new  problem  is  also  generated  using 
information  fi'om  the  previous  discrete  optimization  problem’s  final  solution.  The  inner  and  outer  loop 
structure  defined  for  GHC  algorithms  can  be  used  in  SGHC  algorithms,  where  SGHC  algorithms  restrict 
possible  movement  between  discrete  optimization  problems  to  the  first  iteration  of  the  outer  loop  iterations. 
This  constraint  ensures  that  a  GHC  algorithm  is  applied  to  each  discrete  optimization  problem  at  least  N(k) 
iterations  each  time  it  is  generated  (i.e.,  initially  visited).  Note  that  this  was  not  the  case  for  the 
manufacturing  problem  presented  in  Vaughan  et  al.  (2000),  where  movement  between  discrete  optimization 
problems  was  possible  during  all  iimer  loop  iterations.  The  SGHC  algorithm  is  presented  below  in  pseudo¬ 
code  form. 

Set  the  outer  loop  coimter  boimd  K  and  the  inner  loop  counter  boiuids  N(k),  k  =  0,1,2,. .  .,K 
Define  a  set  of  hill  climbing  (random)  variables  Rt:  Q  x  Q  ->  {-oo,+oo},  k  =  1,2,. .  .,K 

Set  the  iteration  indices  N(0)  =  i  =  0,  k  =  n=  1 

Select  an  initial  discrete  optimization  problem  D(0)  e  S  and  an  initial  solution  (b(0,0)  e  G(0) 

Repeat  while  k  <  K 

Generate  a  discrete  optimization  problem  D(k)  €  Tiset(D(k-l)) 

If  D(k)  ^  D(k-1),  generate  a  solution  ®  e  Q(k)  and  ra(k,  1)  ro 
Else  ffl(k,  1)  <-  ro(k-l,  N(k-1)) 

Repeat  while  n  <  N(k) 

Generate  a  solution  o  g  ri(®(k,  i)) 

Compute  5(ra(k,  i),m)  =  f(ffl)-f(®(k,  i)) 

If  Rk(ro(k,  i),®)  >  5(®(k,  i),®),  then  ®(k,  i+1)  •<-  ® 

If  Rk(®(k,  i),®)  <  5(®(k,  i),®),  then  ®(k,  i+1)  <-  ®(k,  i) 
n  <r-  n+1,  i  i—  i+1 
Until  n  =  N(k) 
n  <-  1,  k  k+1 
Untilk  =  K 

All  SGHC  algorithms  are  formulated  using  three  components,  a  set  of  hill  climbing  random  variables 
{Rk},  a  neighborhood  function  t|  between  solutions,  and  a  neighborhood  function  rjset  between  discrete 
optimization  problems.  The  two-tuple  (k, )  represents  the  inner  loop  iteration  i  =  l,2,...,N(k),  during  outer 
loop  iteration  k  =  1,2,...,K.  D(k)  is  the  discrete  optimization  problem  the  algorithm  is  executing  over 
during  the  k*  outer  loop  iteration,  k  =  1,2,. .  .,K,  where  the  solution  space  of  D(k)  is  depicted  by  Q(k). 

Markov  chain  theory  is  an  effective  tool  for  studying  the  performance  of  local  search  algorithms.  The 
following  analysis  shows  that  an  application  of  the  SGHC  algorithm  can  be  modeled  using  Markov  chains, 
hi  particular,  an  application  of  the  GHC  algorithm  is  first  modeled  using  a  Markov  chain.  Then  an 
application  of  the  SGHC  algorithm  is  modeled  by  a  set  of  Markov  Chains. 

To  show  that  an  application  of  the  GHC  algorithm  can  be  modeled  with  a  Markov  chain,  the  following 
definitions  are  needed.  A  stochastic  process  is  a  family  of  random  variables  defined  on  some  state  space. 
If  there  are  countable  many  members  of  the  family,  the  process  (termed  a  discrete-time  process)  is  denoted 
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by  Qi,  Q2,  where  the  set  of  distinct  values  assumed  by  a  stochastic  process  is  the  state  space.  If  the 
state  space  is  countable  or  finite,  the  process  is  a  chain.  A  stochastic  process  {Qk},  k  =  1,2,...  with  state 
space  Q  =  {©i ,  CO2, . . . }  satisfies  the  Markov  property  if  for  every  n  and  for  all  states  ©  1 ,  ®2>  •  •  •  >  ®n 

Pr{Q„=©„|Q„.i=  ©n-l,  Qn-2=  C0„-2,  ...,Ql=©l}  =  Pr  {Q„=  ©„  |  Qn-1=  ®n-l}  =  Pn(n-1). 

A  discrete-time  stochastic  process  that  satisfies  the  Markov  property  is  a  Markov  chain  (Isaacson  and 
Madsen  1985). 

Let  {Qk}  denote  a  discrete-time  Markov  chain  with  finite  solution  space  Q  =  {©1,  ©2,  ...,  ©pi}.  For 
this  chain  there  are  |Qp  transition  probabilities,  {Py},  i,j  =  1,  2,...,  |Q|.  The  transition  matrix  associated 
with  the  Markov  chain  {Qk}  is  P,  where  Py  is  the  probability  of  moving  from  state  ©i  to  state  ©j. 

An  application  of  a  GHC  algorithm  can  be  modeled  by  a  stochastic  process  {  },  k  =  1,2,...,K,  n  = 

1,2,. .  .,N(k),  e  Q  with  solution  space  Q  =  {©1,  ©2,. . ©pi}  that  satisfies  the  Markov  property  for  every 
n  and  all  states  ©1,  ©2,. . ©n  (i.e.,  {  }  is  a  Markov  chain).  To  see  this,  consider  an  application  of  a  GHC 

algorithm  to  a  discrete  optimization  problem  with  solution  space  Q  =  {©1,  ©2,...,  Wp]}.  Define  gij(k)  to  be 
the  generation  probability  function  for  the  neighborhood  function  p,  where  gy(k)  is  the  probability  that 
©j€r|(©i)  is  generated  during  outer  loop  iteration  k.  Consider  the  inner  loop  iterations  for  fixed  outer  loop 
iteration  k  =  1,  2,...,  K.  Let  { g*  },  k  =  1,  2,...,  K,  n  =  1,  2,...,  N(k),  Ql  e  Q,  be  the  stochastic  process 

where  if  =  ©j,  then  the  GHC  algorithm  is  at  solution  ©j  during  inner  loop  iteration  n  and  outer  loop 

iteration  k  (Johnson  and  Jacobson  2002a,b).  If  the  GHC  algorithm  is  at  solution  ©i  at  inner  loop  iteration  n- 
1,  the  probability  that  the  algorithm  is  at  solution  ©j  at  inner  loop  iteration  n  is 


g,^.(A:)Pr(i?i(©,.,©^.))>5,;,.  for  all  ©,.eQ,  ©^.Gri(©,.),  j*i 


1-  I 

zeJi(a)^) 

z^i 


0 


j  =  i 


otherwise 


independent  of  the  solutions  the  algorithm  visited  at  inner  loop  iterations  1,2,. .  .,n-2.  Therefore 


pr{a'=®ji  e'_2=®i(n-2>-,  ef=®i(i)}=pr{a*=®iia'-.=®i}=Pij(kx 


Hence  the  Markov  property  holds.  Moreover,  for  every  outer  loop  iteration  k,  the  Markov  chain  {Q^}  has 

a  transition  matrix  P(k),  where  Py(k)  is  as  defined  above  . 

Recall,  that  a  SGHC  algorithm  is  applied  to  a  set  of  fundamentally  related  discrete  optimization 
problems.  Movement  between  discrete  optimization  problems  is  only  possible  at  outer  loop  iterations  b=l , 
2,  . . .,  K.  During  the  inner  loop  iterations,  the  SGHC  algorithm  is  executing  over  the  solution  space  of  the 
current  discrete  optimization  problem  using  a  GHC  algorithm 

The  following  analysis  shows  that  for  fixed  outer  loop  iteration  k  =  1,2,...,K,  the  stochastic  process 
corresponding  to  the  SGHC  algorithm  solution  at  iimer  loop  iterations  n  =  1,2,. .  .,N(k)  can  be  modeled  by  a 
Markov  chain  that  corresponds  to  an  application  of  a  GHC  algorithm.  Moreover,  it  is  shown  that  for  outer 
loop  iterations  k  =  1,2,...,K,  the  possible  movement  between  discrete  optimization  problems  is  a  stochastic 
process  that  satisfies  the  Markov  property. 

Consider  an  application  of  a  SGHC  algorithm  to  a  set  of  fundamentally  related  discrete  optimization 
problems  S  =  {Di,  D2,...,  Dm},  where  each  discrete  optimization  problem  Dy,  y  =  l,2,...,m  is  fully  defined 
by  a  solution  space  and  an  objective  function  fy  (i.e.,  Dy=  (Qy,  fy)).  Consider  the  inner  loop  iterations  n 
=  1,2,. . .,N(k),  for  fixed  outer  loop  iteration  k  =  1,2,. .  .,K.  Let  { (Dy)},  k  =  1,2,. . .,K,  n  =  1,2,. . .,N(k)  be 
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the  stochastic  process  where  if  Ql  (Dy)  =  cOj,  then  the  SGHC  algorithm  is  at  solution  ©ke  Qy  at  inner  loop 
iteration  n  of  outer  loop  iteration  k. 

Note  that,  for  all  inner  loop  iterations  n  =  l,2,...,N(k)  of  an  outer  loop  iteration  k  =  1,2,...,K,  the 
SGHC  algorithm  is  executing  over  a  particular  discrete  optimization  problem  from  the  set  of  fundamentally 
related  discrete  optimization  problems  S  =  (Di,  D2,...,  Dm}  using  a  GHC  algorithm.  It  was  shown  that  any 
application  of  a  GHC  algorithm  to  a  discrete  optimization  problem  can  be  modeled  as  a  Markov  chain. 
Therefore,  for  fixed  outer  loop  iteration  k,  the  stochastic  processes  {  Q\  (Dy)},  y  =  1,2,. . .,m,  with  transition 
matrices  Py,  are  the  Markov  chains  that  correspond  to  an  application  of  the  GHC  algorithm  to  the  discrete 
optimization  problems  Dy,  y  =  1,2,. . .,m  for  every  n  =  1,2,. .  .,N(k)  and  for  all  states  ©1,  ©2,. . ©|ny|. 

Movement  between  discrete  optimization  problems  is  a  stochastic  process  that  satisfies  the  Markov 
property.  To  see  this,  define  {'P(k)},  ^(k)  €  S,  k  =  1,2,...  to  be  the  stochastic  process  where  if 'F(k)  =  y, 
then  during  outer  loop  iteration  k,  for  all  inner  loop  iterations  n  =  l,2,...,N(k)  the  SGHC  algorithm  is 
executing  over  solutions  contained  in  the  solution  space  of  the  discrete  optimization  problem  Dy  =  (Qy,  fy). 
If  the  SGHC  algorithm  is  executing  over  Qy  at  outer  loop  iteration  k-1,  then  the  probability  that  the  SGHC 
algorithm  is  executing  over  Qq  during  outer  loop  iteration  k  is 

f'yq(k)  —  (k)  P(I^y»  Dq))? 

independent  of  the  discrete  optimization  problems  the  SGHC  algorithm  visited  at  outer  loop  iterations 
1,2,. .  .,k-2  and  independent  of  all  preceding  inner  loop  iterations.  Therefore, 

Pr {'T(k)  =  q  I  'P(k-1)  =  y,  'P(k-2)  =  yk.2,  . .  .,4^(1)  =  y,}  =  Pr {'P(k)  =  q  |  'P(k-1)  =  y}  =  Tyq(k) 

and  the  Markov  property  holds.  Moreover,  the  Markov  chain  {'P(k)}  has  transition  matrix  T(k),  where 
Tyq(k)  is  as  defined  above. 

Consider  an  application  of  the  SGHC  algorithm  to  a  set  of  fundamentally  related  discrete  optimization 
problems,  S.  This  section  presents  sufficient  conditions  that  guarantee  that  a  SGHC  algorithm  will  (as  k 
approaches  +00)  be  executing  over  the  solution  space  of  each  discrete  optimization  problem  Dy  e  S  with 
probability  1/  jSj,  where  |S|  is  the  cardinality  of  S.  This  result  implies  that,  as  k  approaches  +00,  each 
discrete  optimization  problem  in  S  =  {Di,  D2,. . .,  Dm}  is  being  explored  with  equal  probability. 

This  section  develops  sufficient  conditions  that  place  restrictions  on  the  selection  of  the  problem 
generation  probability  function  (k,  p(Dy,  D,)).  A  discrete  time  Markov  chain  is  a  stationary  Markov 

chain  if  the  probability  of  going  from  one  state  to  another  state  is  independent  of  the  iteration  at  which  the 
transition  is  being  made  (Isaacson  and  Madsen  1985).  That  is,  let  {Xn}  be  a  stationary  Markov  chain  with 
state  space  S={Di,  D2,...,  Dm}.  Then  for  all  states  Dy  and  Dq,  for  all  k  = -(n-1), -(n-2),..., -1,  0,  1,2, ...  , 

Pr  {X„  =  Dy  I  Xn-i  =  Dq}  =  Pr  {X„+k  =  Dy  |  X  „+ic-i  =  Dq} . 

The  long  run  distribution  (stationary  probability  distribution)  of  a  stationary  Markov  chain  with 
corresponding  transition  matrix  T  is  defined hy  11=  ] ,  Tt >  0,  for  all  i  =  1 ,2, . . . ,m,  where 

m 

71  =  ttT  and  =1.  Equivalently,  the  long  run  distribution  of  a  stationary  Markov  chain  is  defined  by  ti  = 

i=l 

[tc^  •••  ;r„],  where 

jij=  lim  Tij^”\ 

K->+oo 

If  Markov  chain  {^(k)}  is  stationary  and  has  a  uniform  long  run  distribution,  then  as  k  approaches 
infinity  the  SGHC  algorithm  is  executing  over  the  solution  space  of  each  discrete  optimization  problem  in 
S  ={Di,  D2,...,  Dm}  with  probability  1  /  m  =  1  /  |S|.  Theorem  2.2  provides  sufficient  conditions  for  the 
selection  of  the  problem  generation  probability  function  (k,  p(Dy,  Dq))  that  guarantee  that  the  Markov 

chain  {'P(k)}  has  a  uniform  long  run  distribution.  Therefore,  when  the  sufficient  conditions  of  Theorem 
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2.2  hold,  the  SGHC  algorithm  will  (as  k  approaches  +oo)  be  in  discrete  optimization  problem  Dy  g  S,  y  = 
1,2,. .  .,m  with  probability  1  /  |S|. 

To  prove  Theorem  2.1,  the  following  definitions  are  needed.  A  subset,  C,  of  the  state  space,  S,  is 
closed  if  Pij  =  0  for  all  i  g  C,  j  g  C.  A  Markov  chain  is  irreducible  if  there  exists  no  nonempty  closed  set 
other  than  S  itself.  If  S  has  a  proper  closed  subset,  it  is  reducible.  State  coi  is  said  to  have  period  if  P  = 

0  whenever  n  is  not  divisible  by  d  and  d  is  the  greatest  integer  with  this  property.  A  state  with  period  one  is 
said  to  be  aperiodic  (Isaacson  and  Madsen  1985). 


Theorem  2.1:  Consider  an  application  of  the  SGHC  algorithm.  Define  (k,  p(Dy,  Dq))  =  (p(Dy5 

Dq))  ~  D  (P(Py>  ^q))>  every  k  =  1,2,...  .  Consider  the  transition  matrix  T  defined  by  Tyq 

~^D  D  (P(I^y’  Dq)).  If  the  transition  matrix  T  is  irreducible  and  aperiodic,  then  the  Markov  chain  {'P(k)} 

has  a  uniform  long  run  distribution.  Moreover,  the  long  run  distribution  of  {'P(k))  is  n  =  [1/|S|  1/|S|  ... 
1/|S|]. 

Proof:  See  Vaughan  and  Jacobson  (2003b). 

Theorem  2.1  shows  that  if  the  problem  generation  probability  function  is  chosen  such  that  the  Markov 
chain  {^(k)}  is  stationary  and  the  associated  transition  matrix  is  symmetric,  then  the  Markov  chain  {'P(k)} 
has  a  uniform  stationary  distribution.  The  second  set  of  sufficient  conditions.  Theorem  2.2,  allows  for  the 
development  of  a  nonstationary  Markov  chain.  To  present  Theorem  2.2,  a  discussion  of  weak  and  strong 
ergodicity.  Lemmas  2.1  and  Lemma  2.2  are  needed. 

For  finite  nonstationary  Markov  chains,  the  transition  matrices  T(k)  that  contain  the  probabilities  of 
moving  from  state  Dy  to  state  Dq  at  outer  loop  iteration  k,  are  functions  of  k.  To  define  weak  and  strong 
ergodicity  of  nonstationary  Markov  chains,  several  definitions  are  needed.  Define  the  one  norm  of  a  vector 

m  m 

f-(fb  f2v5  fm)  by  II  f  II  =XI  fi  I  •  Define  the  infinity  norm  of  matrix  T(k)  by  ||T(k)||  =max  XI  ^7 1  (Atkinson 

1=1  i  7=1 

1989).  Let  T(l),  T(2),  ...  be  the  transition  matrices  for  a  nonstationary  Markov  chain  with  starting  vector 
Define  f^’'‘^=f<°^0'+l)T(j+2) ...  T(k). 

A  nonstationary  Markov  chain  is  weakly  ergodic  if,  for  all  j, 

lim  sup  ||f^'’">-g«’''>||  =  0, 

«->+«  y'(0)g(0) 

where  f^°^  and  g^°^  are  starting  vectors.  A  nonstationary  Markov  chain  is  strongly  ergodic  if  there  exists  a 
vector  q  =  (qj,  q2, . . qm),  with  ||  q  ||=1  and  qi>0,  for  i=l,  2, . . .,  m  such  that,  for  all  j, 

lim  sup  =  0, 

M-^+QO  y.(0) 


where  is  a  starting  vector  (Isaacson  and  Madsen  1985).  The  following  result  from  Isaacson  and  Madsen 
(1985)  is  needed  in  the  proof  of  Lemma  2.2 

Lemma  2.1:  Let  {T(k)}  be  a  sequence  of  transition  matrices  corresponding  to  a  nonstationary  weakly 
ergodic  Markov  chain.  If  there  exists  a  corresponding  sequence  of  left  eigenvectors  {7i(k)},  for  {T(k)}, 
satisfying 

+00 

ZII’t(*)-<*+i)ll<'' 

k=\ 

then  the  chain  is  strongly  ergodic  and  for  every), 

lim  sup  llf^’"^-  Till  =  0, 

«->+00  y.(0) 

where 


lim  %  (k)  =  n. 
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Proof:  See  Isaacson  and  Madsen  (1985). 

Lemma  2.2  shows  that  if  the  problem  generation  probability  function  is  selected  such  that  the 
corresponding  nonstationary  Markov  chain  {'P(k)},  'P(k)  e  S,  k  =1,2,...  is  weakly  ergodic  and  the 
transition  matrices  T(k)  are  S3nnmetric  for  every  k,  then  the  nonstationary  Markov  chain  is  strongly  ergodic. 


Lemma  2.2:  Consider  an  application  of  the  SGHC  algorithm.  Assume  that  the  nonstationary  Markov  chain 
{'F(k)},  'P(k)  e  S,  k  =  1,2,...  is  weakly  ergodic  and  that  hj^  ^  (k,  p(Dy,  Dq))  =  ^  (k,  p(Dq,  Dy)),  for 


every  k.  Then  the  nonstationary  Markov  chain  {'P(k)}  is  strongly  ergodic  and  for  every), 

lim  sup  ||f’’"^-n||  =  0, 

«^+CO  y.(0) 


where  7t=[l/|S|  1/|S| ...  1/|S|]. 

Proof;  See  Vaughan  and  Jacobson  (2003b). 


When  the  assumptions  in  Lemma  2.2  hold,  then  the  SGHC  algorithm  will  (as  k  approaches  +oo)  be 
executing  over  the  solution  space  of  each  discrete  optimization  problem  contained  in  the  set  of 
fundamentally  related  discrete  optimization  problems  with  equal  probability.  Theorem  2.2  shows  that  this 
result  implies  that  for  every  s  >  0,  there  exists  an  outer  loop  iteration  such  that  for  all  future  outer  loop 
iterations,  the  SGHC  algorithm  is  executing  over  the  solution  space  of  each  discrete  optimization  problem 
in  the  set  S  with  probability  1/|S|  ±  e. 


Theorem  2.2:  Consider  an  application  of  the  SGHC  algorithm.  Assume  that  the  Markov  chain  {'P(k)}, 
'P(k)  G  S,  k  =  1,2,...  is  weakly  ergodic  and  that  hj^  ^  (k,  p(Dy,  Dq))  =  (k,  p(Dq,  Dy)),  for  every  k  = 

1,2,. . .  and  for  every  q,y  =  1,2,. .  .,m  (i.e.,  the  corresponding  transition  matrix  is  symmetric).  Then  for  every 
8  >  0,  there  exists  a  k(8)  e  7^,  such  that  1/|S|  -  8  <  Pr{T'(  k)  =  y}  <  1/|S|  +  s,  for  all  outer  loop  iterations  k  > 
k(8)  and  for  every  Dy  g  S,  y  =  1,2,. .  .,m. 

Proof:  See  Vaughan  and  Jacobson  (2003b). 

Theorem  2.2  establishes  that  the  SGHC  algorithm  can  be  designed  to  ensure  that  the  probability  of 
visiting  each  problem  is  uniformly  distributed.  The  following  discussion  illustrates  this  result  on  a  set  of 
fundamentally  related  discrete  optimization  problems  using  the  traveling  salesman  problem.  The  traveling 
salesman  problem  (TSP)  is  a  well-known  NP-hard  discrete  optimization  problem  (Lawler  et  al.  1985).  The 
TSP  is  used  to  illustrate  various  local  search  algorithms  because  it  is  useful  for  modeling  a  variety  of  real 
world  problems.  For  example,  traditional  applications  of  the  TSP  include  a  variety  of  vehicle  routing  and 
scheduling  problems.  More  recently,  applications  of  the  TSP  have  been  expanded  to  include  modem 
applications  like  the  printing  of  circuit  boards,  x-ray  crystallography,  overhauling  of  gas  turbine  engines, 
and  the  controlling  of  industrial  robots  (Johnson  and  Jacobson  2002a, b). 

To  formally  describe  the  TSP  (Lawler  et  al.  1985),  define  a  graph  to  be  a  finite  set  of  vertices,  some 
pairs  of  which  are  joined  by  edges.  A  cycle  in  a  graph  is  a  set  of  vertices  of  the  graph,  which  is  such  that  it 
is  possible  to  move  from  vertex  to  vertex,  along  edges  of  the  graph,  so  that  all  vertices  are  encountered 
exactly  once,  finishing  at  the  start.  If  a  cycle  contains  all  the  vertices  of  the  graph,  it  is  called  a 
Hamiltonian  cycle  (or  tour).  The  TSP  is  defined  as  follows  (Garey  and  Johnson  1979). 


Traveling  Salesman  Problem  (TSP) 

Instance:  Given  a  set  of  n  cities  C  =  {ci,  C2,...,  Cn}  and  a  distance  matrix  D  that  represents  the  cost  of 
traveling  between  the  cities  in  the  set  C. 

/J-l 

Question:  Find  a  Hamiltonian  cycle  H  =  (c’i,c’2,. .  .,c’n)  that  minimizes  f(H)  =  X  f^(c’j,c’j+i)  +  D(c’„,  c’l). 

y=i 

An  instance  of  a  TSP  is  a  discrete  optimization  problem,  where  the  solution  space  is  the  set  of  possible  all 
Hamiltonian  cycles  (with  each  tour  consisting  of  n  cities),  Q  =  {coj,  (i>2,  C0(n.i)!/2} •  The  objective  function 

value  for  each  solution  C0i=(c’i,  c’2,  ...,  c’n)  e  G  is  the  sum  of  the  distances  the  tour  depicts,  f(c0i)  = 
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n-\ 

2  -D(c’j,  c’j+i)  +  D(c’n,  c’l).  The  optimal  objective  function  value  represents  the  shortest  distance 
y=i 

traveled. 

The  Multiple  Traveling  Salesman  Problem  (MTSP)  is  introduced  to  illustrate  the  application  of  SGHC 
algorithms.  The  MTSP  is  defined  as  follows. 

Multiple  Traveling  Salesman  Problem  (MTSP) 

Instance:  Given  a  set  of  n  cities  Ob  =  {ci,  Ci,---,  Cn},  a  set  of  m  subsets  of  Ob,  O  =  {Ci,  Ci,--.,  Cm},  and  a 
distance  matrix  D  that  represents  the  cost  of  traveling  between  the  cities  in  the  set  Ob. 

Question:  Find  a  Hamiltonian  cycle  H  =  (c’l,  c’2,...,  where  there  exists  a  Cy,  y  =  l,2,...,m  such 

n~\ 

that  c’j  €  Cy,  for  every  j  =  1,2,...,  n,  and  f(H)  =  Y.  -DCc’j,  c’j+j)  +  c’l)  is  minimized. 

M 

Note  that  each  of  the  sets  Cy  e  O  represents  an  instance  of  the  TSP,  Dy.  Since  the  TSP  can  be 
formulated  as  a  discrete  optimization  problem,  then  the  MTSP  problem  can  be  represented  by  set  of 
discrete  optimization  problems  S  =  {Di,  D2, ...,  Dm}.  The  set  S  =  {Di,  D2, ...,  Dm}  is  a  set  of  fundamentally 
related  discrete  optimization  problems.  To  see  this,  note  that  each  discrete  optimization  problem  Dy  e  S,  y 
=  1,2,. ..,m  is  fully  defined  by  CyC  06  =  {ci,  C2,  ...,  Cn}.  Therefore,  Cy  is  the  fundamental  relation  set  of 
discrete  optimization  problem  Dy. 

The  following  computational  results  illustrate  how  a  MTSP  can  be  addressed  with  SGHC  algorithms. 
For  comparison  purposes,  GHC  algorithms  are  also  applied  to  the  individual  members  in  the  set  of 
fundamentally  related  discrete  optimization  problems.  These  computational  results  suggest  that 
optimal/near  optimal  solutions  can  be  reached  in  less  total  iterations  using  a  SGHC  algorithm. 

To  develop  a  MTSP,  a  set  consisting  of  20  cities  was  generated  by  randomly  locating  each  city  on  a 
1000  by  1000  unit  grid.  Four  TSPs  of  size  18  were  generated  by  randomly  selecting  18  of  the  20  cities  for 
each  traveling  salesman  problem.  In  the  case  where  an  identical  set  of  cities  was  generated  for  two  or  more 
of  the  problems,  a  completely  new  set  of  discrete  optimization  problems  was  generated,  resulting  in  four 
distinct  randomly  generated  TSPs. 

The  four  randomly  generated  TSPs  are  arbitrarily  labeled  Di,  D2,  D3,  and  D4.  The  detachment  metric 
p(Dy,  Dq),  between  the  TSPs,  was  calculated  for  all  y,q  =  1,2, 3, 4.  The  distance  matrix  and  distance 
diagram  are  depicted  in  Figure  2.1. 

Figure  2.1:  Distauce  Matrix  and  Distance  Diagram 
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Computational  results  with  Monte  Carlo  search,  pure  local  search,  and  simulated  aimealing  using 
SGHC  algorithms  are  reported.  For  comparison  purposes,  computational  results  with  Monte  Carlo  search, 
pure  local  search,  and  simulated  annealing  using  GHC  algorithms  are  also  reported.  The  2-Opt 
neighborhood  function  was  used  for  all  executions  of  the  SGHC  and  GHC  algorithms.  For  the  SGHC 
algorithms,  the  problem  generation  probability  function  is  defined  as 

6^  n  (k,  p(Dy,  Dq))  =  [l/p(Dy,  Dq)]  /  [  t  (l/p(Dy,  D;))],  y  q 

^  /=! 

i^y 
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and 


^  (k,  p(Dy,  Dy))  =  1-  E  [l/(p(Dy,  Dq)]  /  [  i  (l/p(Dy,  Di))]  =  0, 

^  ^  «=1  (=1 
q*y  i*y 

for  every  y,q  =  1,2,3 ,4,  y  q  for  every  k  =  1,2,. .  .,K. 

This  problem  generation  probability  function  is  defined  such  that  the  associated  transition  matrix  is 
symmetric  and  the  Markov  chain  {'F(k)}  is  stationary.  Therefore,  by  Theorems  2.1  and  2.2,  the  Markov 
chain  {'F(k)}  has  a  uniform  stationary  distribution,  hence  as  k  approaches  infinity,  the  SGHC  algorithm  is 
executing  over  the  solution  space  of  each  discrete  optimization  problem  in  S  =  {Di,  D2,...,  Dm}  with 
probability  1/m  =  1/|S|  =  1/4.  Moreover,  this  problem  generation  function  guarantees  the  discrete 
optimization  problem  over  which  the  SGHC  algorithm  is  executing  changes  at  every  outer  loop  iteration  k 
(i.e.,  TO  ^  'F(k-l),  for  all  k  =  1,2,. . .). 

Executions  with  different  values  of  K  and  N  =  N(k),  k  =  1,2,... ,K  are  reported.  To  compare  the 
performance  of  applying  a  SGHC  algorithm  versus  applying  a  GHC  algorithm,  the  inner  and  outer  loop 
bounds  of  the  SGHC  algorithm  were  doubled.  Therefore,  the  total  number  of  iterations  that  the  SGHC 
algorithm  executes  is  equal  to  the  total  number  of  iterations  executed  using  the  GHC  algorithm  for  the  four 
individual  problems.  Let  R  e  represent  the  total  number  of  replications  executed  for  each  SGHC  and 
GHC  algorithm  formulation.  For  each  replication,  a  different  randomly  generated  initial  tour  was  used. 
The  means,  p,  standard  deviations,  ct,  and  the  minimum  and  maximum  distances,  were  computed  from  the 
optimal  tour  distances  across  these  R  replications.  The  value  y  in  Tables  2.1  to  2.6  represents  the  number 
of  replications  for  which  the  algorithms  find  the  minimum  distance  tour.  For  simulated  annealing,  tk  is 
updated  by  multiplying  the  previous  temperature  parameter  by  the  increment  multiplier  p  =  .90  (i.e.,  tk  = 
Ptk-i),  with  initial  temperature  parameter  to  =  2000. 


Table  2.1;  GHC  Algorithm  Results;  Pure  Local  Search 


Outer  and  Inner  Loop  Bounds 

y/R 

a 

Minimum 

Maximum 

K=100,  M=100 

4/30 

3864.3 

36.9726 

3805.4 

3916.0 

K=200,  M=100 

3863.4 

32.7691 

3907.3 

K=300,  M=75 

3876.6 

36.9549 

3953.7 

K=400,  M=75 

3885.3 

42.1893 

3973.4 

K=400,  M-50 

3867.9 

37.7520 

3916.0 

K=800,  M=50 

1/15 

3872.8 

35.2469 

3805.4 

3916.0 

Table  2.2;  SGHC  A1 

gorithm  Results;  Pure  Local  Search 

Outer  and  Inner  Loop  Bounds 

y/R 

a 

Minimum 

Maximum 

K=100,  M=100 

10/30 

3819.4 

13.2943 

3805.4 

3831.8 

K=200,  M=100 

9.0535 

3831.6 

K=300,  M=75 

9.0535 

3831.6 

K=400,  M=75 

24/30 

3807.2 

6.6434 

3831.6 

K=400,  M=50 

9/30 

3818.5 

13.3165 

3831.6 

K=800,  M=50 

12/15 

3810.7 

10.8417 

3805.4 

3831.6 

The  results  in  Table  2.1  suggest  that  when  the  number  of  outer  loop  iterations  for  a  pure  local  search 
GHC  algorithm  is  increased  from  100  to  200,  performance  of  the  algorithm  shows  no  improvement. 
However,  Table  2.2  suggests  that  the  performance  of  a  pure  local  search  SGHC  algorithm  improves 
significantly  (as  measured  by  p)  when  the  number  of  outer  loop  iterations  is  increased  from  100  to  200. 
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Table  2.3;  GHC  Algorithm  Results;  Simulated  Annealing 


Outer  and  Inner  Loop  Bounds 

to 

y/R 

11 

a 

Minimum 

Maximum 

IBtM 

41.4192 

3805.4 

3953.7 

35.7662 

3805.4 

3916.0 

K=300,  M=75 

K=400,  M=75 

UtMJI 

3805.4 

3907.5 

K=400,  M=50 

2000 

3/30 

3868.8 

39.4863 

3805.4 

3916.0 

K=800,  M-50 

2000 

1/15 

3885.9 

29.7020 

3805.4 

3916.0 

Table  2.4;  SGHC  Algorithm  Results;  Simulated  Annealing 


Outer  and  Inner  Loop  Boimds 

To 

y/R 

a 

Minimum 

Maximum 

K=100,  M=100 

16/30 

3814.1 

12.5549 

3805.4 

3831.6 

K=200,  M=100 

19/30 

3808.0 

7.9899 

3805.4 

3831.6 

K=300,  M=75 

20/30 

3808.9 

9.0535 

3805.4 

K=400,  M=75 

23/30 

3808.9 

9.0535 

3805.4 

K=400,  M=50 

2000 

7/30 

3831.6 

13.1248 

3805.4 

3831.6 

K=800,  M=50 

2000 

1/15 

3813.6 

12.5352 

3805.4 

3831.6 

Table  2.5: 

GHC  Algorithm  Resul 

ts;  Monte  Carlo  Search 

Outer  and  Inner  Loop  Bounds 

y/R 

a 

Minimum 

Maximum 

K=100,  M=100 

6390.0 

294.6998 

5634.9 

6805.9 

K=200,  M=100 

6195.6 

239.7325 

5575.7 

6656.3 

K=300,  M=75 

mi 

5293.0 

6613.9 

K=400,  M=75 

1/30 

292.5790 

5310.2 

6488.3 

K=400,  M=50 

1/30 

6122.6 

287.3693 

5291.7 

6575.7 

K=800,  M=50 

1/15 

6007.1 

231.746 

5479.1 

6328.1 

Table  2.6;  SGHC  Algorithm  Results;  Monte  Carlo  Search 


Outer  and  Inner  Loop  Bounds 

y/R 

a 

Minimum 

Maximum 

K=100,  M=100 

6758.2 

202.9184 

5793.3 

6758.2 

K=200,  M=100 

6109.3 

301.5243 

5243.7 

6611.5 

K=300,  M=75 

1/30 

6214.6 

204.3092 

K=400,  M=75 

1/30 

6054.3 

5614.9 

6462.8 

K=400,  M=50 

1/30 

6265.8 

5877.3 

6659.3 

K=800,  M=50 

1/15 

5954.7 

301.7151 

5299.6 

6382.3 

The  results  in  Tables  2.5  and  2.6  suggest  that  there  is  little  difference  in  performance  between  Monte 
Carlo  search  GHC  algorithms  and  Monte  Carlo  search  SGHC  algorithms.  Note  that  this  result  occurs  since 
Monte  Carlo  search  accepts  neighboring  solutions  independent  of  their  objective  function  value. 
Therefore,  when  movement  between  discrete  optimization  problems  occurs  there  is  no  exchange  of 
common  information  between  the  discrete  optimization  problems.  Overall,  the  results  in  Tables  2.1 
through  2.4  suggest  that  the  SGHC  algorithms  outperform  the  GHC  algorithms.  The  minimum  distance 
found  over  the  R  replications  using  SGHC  algorithms  is  significantly  smaller  than  the  minimum  distance 
found  over  the  R  replications  using  GHC  algorithms  for  both  the  simulated  annealing  and  pure  local  search 
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algorithms.  Additionally,  the  standard  deviation  of  the  optimal  values  over  the  R  replications  is  smaller 
using  the  SGHC  algorithms. 

Figures  2.2  through  2.4  depict  plots  comparing  the  performance  of  the  SGHC  algorithms  and  the  GHC 
algorithms.  To  obtain  this  data,  fifteen  replications  of  each  SGHC  algorithm  and  GHC  algorithm 
formulation  were  executed.  For  each  replication,  a  different  randomly  generated  initial  solution  was  used. 
The  mean  of  the  optimal  distances  across  the  fifteen  replications  for  the  GHC  algorithm  are  plotted  with  a 
solid  blue  line.  The  standard  deviations  of  the  optimal  distances  for  the  GHC  algorithm  across  the  fifteen 
replications  are  plotted  with  a  dashed  blue  line.  The  means  of  the  optimal  distances  across  the  fifteen 
replications  for  the  SGHC  algorithm  are  plotted  with  a  solid  red  line.  The  standard  deviations  of  the 
optimal  distances  for  the  SGHC  algorithm  across  the  fifteen  replications  are  plotted  with  a  dashed  red  line. 
The  number  of  outer  loop  iterations  executed  was  800  and  the  number  of  inner  loop  iterations  executed  was 
50  for  every  formulation.  For  simulated  annealing,  tk  is  updated  by  multiplying  the  previous  temperature 
parameter  by  the  increment  multiplier  P=.90  with  initial  temperature  parameter  to==2,000. 

Figure  2.2:  Pure  Local  Search 

Pure  Local  Search 


Figure  2.3:  Simulated  Annealing 

Simulated  Annealing 


GHC 

GHC  stdev 
SGHC 
SGHC  stdev 


GHC 

GHC  stdev 
SGHC 
SGHC  stdev 


16 


Figure  2.4:  Monte  Carlo  Search 


Monte  Carlo  Search 
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Figures  2.2  and  2.3  suggest  that  the  SGHC  algorithms  outperform  the  GHC  algorithms,  as  measured  by 
|Li.  The  minimum  distance  found  over  the  fifteen  replications  using  SGHC  algorithms  is  significantly 
smaller  after  2500  iterations  than  the  minimum  distance  found  over  the  fifteen  replications  using  GHC 
algorithms  for  both  the  simulated  annealing  and  pure  local  search  algorithms.  Since  the  SGHC  algorithms 
periodically  move  between  discrete  optimization  problems  they  do  not  get  trapped  at  local  minima  as 
frequently  as  GHC  algorithms.  Moreover,  the  standard  deviation  band  of  the  optimal  values  over  the 
fifteen  replications  is  smaller  using  SGHC  algorithms.  Figure  2.4  suggests  that  there  is  no  significant 
difference  in  the  performance  of  Monte  Carlo  Search  SGHC  and  GHC  algorithms.  As  noted  previously, 
this  result  occurs  since  Monte  Carlo  search  accepts  neighboring  solutions  independent  of  their  objective 
function  value.  Therefore,  when  movement  between  discrete  optimization  problems  occurs  there  is  no 
exchange  of  common  information  between  the  discrete  optimization  problems. 

Figures  2.5  through  2.7  depict  plots  comparing  the  performance  of  the  SGHC  algorithms  and  the  GHC 
algorithms.  To  obtain  this  data,  thirty  replications  of  each  SGHC  algorithm  and  GHC  algorithm 
formulation  were  executed.  For  each  replication,  a  different  randomly  generated  initial  solution  was  used. 
The  mean  of  the  optimal  distances  across  the  thirty  replications  for  the  GHC  algorithm  are  plotted  with  a 
solid  blue  line.  The  standard  deviations  of  the  optimal  distances  for  the  GHC  algorithm  across  the  thirty 
replications  are  plotted  with  a  dashed  blue  line.  The  means  of  the  optimal  distances  across  the  thirty 
replications  for  the  SGHC  algorithm  are  plotted  with  a  solid  red  line.  The  standard  deviations  of  the 
optimal  distances  for  the  SGHC  algorithm  across  the  thirty  replications  are  plotted  with  a  dashed  red  line. 
The  number  of  outer  loop  iterations  executed  was  400  and  the  number  of  inner  loop  iterations  executed  was 
75  for  every  formulation.  For  simulated  annealing,  4  is  updated  by  multiplying  the  previous  temperature 
parameter  by  the  increment  multiplier  p  =  .90  with  initial  temperature  parameter  to  =  2,000. 


17 


Figure  2.5:  Pure  Local  Search 
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Figure  2.6:  Simulated  Annealing 
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Figure  2.7:  Monte  Carlo  Search 
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Figures  2.5  and  2.6  suggest  that  the  SGHC  algorithms  outperform  the  GHC  algorithms,  as  measured  by 
p.  The  minimum  distance  found  over  the  thirty  replications  using  SGHC  algorithms  is  significantly  smaller 
after  2500  iterations  than  the  minimum  distance  found  over  the  thirty  replications  using  GHC  algorithms  for 
both  the  simulated  annealing  and  pure  local  search  algorithms.  Since  the  SGHC  algorithms  periodically 
move  between  discrete  optimization  problems  they  do  not  get  trapped  at  local  minima  as  frequently  as  GHC 
algorithms.  Moreover,  the  standard  deviation  band  of  the  optimal  values  over  the  thirty  replications  is 
smaller  using  SGHC  algorithms.  Figure  2.7  suggests  that  there  is  no  significant  difference  in  the 
performance  of  Monte  Carlo  search  SGHC  and  GHC  algorithms.  As  before,  this  result  occurs  since  Monte 
Carlo  search  accepts  neighboring  solutions  independent  of  their  objective  function  value.  Therefore,  when 
movement  between  discrete  optimization  problems  occurs  there  is  no  exchange  of  common  information 
between  the  discrete  optimization  problems. 

The  SGHC  algorithm  is  a  new  approach  for  addressing  a  set  of  fundamentally  related  discrete 
optimization  problems  that  can  be  more  efficient  than  the  traditional  approach  of  addressing  each  discrete 
optimization  problem  in  the  set  S  individually  with  a  local  search  algorithm.  For  example,  SGHC 
algorithms  allow  practitioners  to  make  a  single  algorithm  run  over  a  set  of  fundamentally  related  discrete 
optimization  problems.  Moreover,  the  computational  results  presented  suggest  that  a  SGHC  algorithm  can 
outperform  the  GHC  algorithm.  The  development  of  the  SGHC  algorithm  and  the  mathematical  results  in 
this  paper  make  it  possible  for  the  SGHC  algorithm  to  be  adapted  and  used  to  approach  a  variety  of  real- 
world  problems. 


3.  P-Acceptable  Solution  Performance  Measures 

An  important  breakthrough  has  been  the  introduction  and  development  of  (3-acceptable  solution  analysis 
for  GHC  algorithms.  To  obtain  these  new  results,  several  important  research  development  milestones  had 
to  be  attained.  All  these  results  are  reported  in  Orosz  and  Jacobson  (2002a,b). 

The  current  literature  on  asymptotic  convergence  properties  and  finite-time  performance  measures 
focuses  primarily  on  reaching  a  global  minimum.  However,  in  practice,  solutions  that  are  close  enough 
to  a  global  minimum  are  often  acceptable.  Define  solutions  that  have  objective  function  value  no  greater 
than  P  as  P-acceptable  solutions,  where  p  denotes  the  maximum  acceptable  objective  function  value 


19 


(necessarily  larger  than  the  global  minimum  objective  function  value).  This  paper  analyzes  the  finite¬ 
time  behavior  of  GHC  algorithms  in  reaching  P-acceptable  solutions. 

To  illustrate  how  the  P-acceptable  solution  framework  can  be  applied,  a  variation  of  the  SA 
algorithm  that  uses  a  fixed  cooling  schedule  term  static  simulated  annealing  (S^A)  is  used.  Desai  (1999) 
and  Cohn  and  Fielding  (1999)  present  results  regarding  the  finite-time  behavior  of  the  S^A  algorithm  that 
suggest  the  probability  of  visiting  a  globally  optimal  solution  is  maximized  by  keeping  the  temperature 
fixed  rather  than  allowing  the  temperature  to  approach  zero.  Fielding  (2000)  presents  an  insightful 
empirical  study  that  suggests  that  there  is  an  optimal  fixed  temperature  for  S^A  for  different  problems.  A 
disadvantage  of  maintaining  a  fixed  temperature  is  that  the  necessary  convergence  condition  (the  cooling 
schedule  must  approach  zero  as  the  number  of  iterations  approaches  infinity,  e.g.,  see  Hajek  1988)  for 
SA,  is  violated.  During  an  execution  of  S^A,  the  algorithm  may  be  in  a  neighborhood  of  a  global 
optimum,  yet  not  visit  this  global  optimum  since  there  exists  a  positive  probability  of  not  visiting  it  (i.e., 
hill  climbing).  An  advantage  of  S^A  is  that  there  is  a  nondecreasing  positive  probability  of  escaping  a 
local  optimum  throughout  the  algorithm’s  execution. 

This  research  looks  at  finite-time  performance  measures  for  identifying  y9-acceptable  solutions.  In 
particular,  the  paper  derives  expression  for  the  expected  number  of  iterations  to  visit  the  set  of  p- 
acceptable  solutions,  and  uses  these  expressions  to  obtain  upper  and  lower  bounds  for  this  expectation 
that  can  be  estimated  using  finite  length  runs.  These  bound  estimators  are  then  computed  from  empirical 
data  generated  from  independent  replications  of  (moderately  short)  S^A  algorithm  runs.  The  resulting 
estimators  and  the  associated  estimation  procedure  provide  information  that  can  be  used  by  practitioners 
to  evaluate  how  long  an  S^A  algorithm  should  be  run  (as  measured  by  the  expected  number  of  iterations) 
to  be  able  to  visit  the  set  of  y9-acceptable  solutions  for  different  values  of  p.  Note  that  the  estimation 
procedure  and  the  need  for  multiple  replications,  as  in  its  current  form,  suggest  that  it  may  not  be 
practical  to  be  used  for  very  large  problem  instances. 

The  objective  function  and  the  neighborhood  function  for  a  discrete  optimization  problem  allows  the 
solution  space,  Q,  to  be  decomposed  into  two  mutually  exclusive  and  exhaustive  sets: 

-  the  set  of  globally  optimal  solutions,  G  e  ^ficd)  for  all  g  Q}, 

-  the  set  of  all  other  solutions,  G'  =  {tu  g  Q:y(ty*)  <j{cd),co*  e  G}  where  G  u  G‘'=  Q. 

In  many  applications,  the  goal  when  addressing  a  discrete  optimization  problem  is  to  identify  a  globally 
optimal  solution  a*  g  G.  However,  from  a  practical  point  of  view,  solutions  that  are  close  enough  to  a 
globally  optimal  solution  (where  close  enough  is  measured  in  terms  of  the  objective  function  value)  for  a 
discrete  optimization  problem  may  be  acceptable.  To  describe  such  solutions,  define  the  set  of  P- 
acceptable  solutions 

{coeQ.\j{a)<p},pe  9i..  (1) 

Note  that  if  P  <fico*),  m*  g  G,  then  Dp  =  0.  Moreover,  \f  P  >  max^  e  n  J{cd),  then  Dp  =  Cl.  Lastly, 

lim  Dp  =  G,  hence  G  is  the  upper  (right)  limit  of  Dp  as  p  approaches^^y*)  from  above. 
p^f(<oy 

Simulated  annealing  (SA)  is  a  particular  GHC  algorithm,  with  hill  climbing  random  variable, 
Rii(o,a)’)  =  -T(k)*ln(l-Uk),  co’eTj(cd),  where  the  {Uk}  are  HD  U(0,1)  random  variables  and  T(k), 
k=l,2,...,K,  are  the  temperate  parameters  that  define  the  cooling  schedule.  Therefore,  -f{o)(i))  > 

0,  then  cop)  is  accepted  as  the  current  solution  with  probability  The  temperature 

parameters  are  generally  set  such  that  they  gradually  decrease  to  zero.  This  means  that  as  k  approaches 
infinity,  the  SA  algorithm  behaves  similar  to  pure  local  search  algorithm,  hence  it  will  eventually 
terminate  at  either  a  global  or  a  local  optimum. 

Static  simulated  annealing  (S^A)  is  a  particular  type  of  SA  algorithm,  where  the  temperature  is  held 
fixed  (or  static)  during  the  entire  algorithm  execution.  If fico)  -fpop))  >  0,  then  cop)  is  accepted  as  the 
current  solution  with  probability  ^  where  T  =  T(k),  k  =  1,2,... ,K  is  held  constant  at  each 
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iteration.  The  S^A  and  SA  pseudo-code  are  identical,  except  that  only  a  single  temperature  T  is  required 
for  S^A  throughout  the  entire  algorithm  execution. 

The  one-step  p-acceptable  probability  is  now  introduced  and  described  as  a  finite-time  performance 
measure  for  GHC  algorithms.  Define  A  to  be  a  GHC  algorithm  applied  to  an  instance  of  a  discrete 

K 

optimization  problem,  where  h  =  X  N(k)  represents  the  total  number  of  algorithm  iterations  executed. 

i=l 

Assume  that  for  algorithm  A,  Rk{(o(i),co)  >  0,  o>(i)  €  fl,  m  e  t]{co(i)),  for  all  iterations  i  =  1,2,. .  .,h,  and  for 
all  outer  loop  iterations  k  =  1,2,...,K.  At  each  iteration  h,  define  the  event  D(h,yS)  on  the  probability 
space  (S,  3,  P),  termed  the  set  of  P-acceptable  solutions,  as 

DA(h,yff)  =  D(h,)ff)  ={{co(l),a)(2),....co(h)y.  co(i)  e  i=  l,2,...,h,fio)(i))  <p  for  some  f  =  l,2,...,h} 
={{co(l),co(2),...,co(h)y.  co(i)  g  Q.,i=  l,2,...,h,  tuflj  g  Dyjforsomei=  1,2,. ..,h},  (2) 

where  p  is  an  objective  function  value  threshold.  Therefore,  the  complementary  event  is 

DA‘'(h,y9)  =  D‘'(h,;9)  =  {(co(l),co(2),...,(o(h)):  co(i)  eQ,i=  2,..., h, f{co(i))  > P  for  a\\  i=  l,2,...,h} 

=  {(co(l),co(2),...,co(h)):  co(i)  e  Q,i=  1,2,. ..,h,  co(i)  g  D^sfor  all  i=  1,2,. ..,h}.  (3) 

The  event  D(h,y9)  defines  sequences  of  h  solutions  that  result  from  the  execution  of  algorithm  A  over 
h  iterations,  where  one  or  more  solutions  have  objective  function  values  less  than  or  equal  to  p.  The 
definition  of  D(h,y9)  implies  that  D(i,;^  c  D(i+l,y5),  for  all  iterations  i  =  l,2,...,h,  hence  {D(h,^)}  is  a 
telescoping,  non-decreasing  sequence  of  events  in  h.  From  (2)  and  (3),  the  one-step  p-acceptable 
probability  is  defined  as 

r(j,y9)  =  P{D(j,y9)|D‘=0-l,y9)}  =V{{co(l),co(2),...,(D(j)y.  co{i)  g  Q.,  i=  1,2,... fico(i))> p^ox^]M  = 

\,2,...,y\,Aco(]))  <P}IV{m-hP)}  (4) 

The  one-step  yS-acceptable  probability  at  iteration  j  provides  a  finite-time  performance  measure  for 
the  effectiveness  of  a  GHC  algorithm,  namely  the  ability  of  the  algorithm  to  visit  an  element  of  the 
solution  space  with  objective  function  value  less  than  or  equal  to  p  at  iteration  j  given  that  it  has  not 
already  visited  such  a  solution  over  the  first  j-1  iterations. 

The  one-step  y9-acceptable  probability  can  be  used  to  obtain  a  closed  form  expression  for  (2).  Lemma 
3.1  captures  this  relationship. 

Lemma  3.1:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that  P{D°(0,y8)}  = 

1.  ThenP{DM}  =  l-fl  [l-r(j,A)] 

Proof.  See  Orosz  and  Jacobson  (2002a). 

The  one-step  y9-acceptable  probability  can  be  used  to  express  the  expected  number  of  iterations  to  reach 
the  ^-acceptable  set  for  the  first  time.  In  particular,  define  the  random  variable  to  represent  the 
minimum  number  of  iterations  needed  to  reach  an  element  in  the  set  of  y9-acceptable  solutions, 

=  min  (j  >  1  :  f(coQ))  <  P) .  (5) 

The  relationship  between  and  the  one-step  ^-acceptable  probability  is  described  in  Lemma  3.2. 

Lemma  3.2:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that 
P{D"(0,y9)}=l.  Then 

(a)  P{x^  >  h}  =n  [1  - r(jM  =  P{D'(h,  P)} 

M 

+00 

(b)  P{t^  <  +00}  =  1  -  n  [1  -  T(i,p)] 
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(c)  P{x^  =  h}=r(li,yS)n  [l-rQM- 

7=1 

Proof:  See  Orosz  and  Jacobson  (2002a). 

Theorem  3.1  provides  an  expression  for  the  conditional  expectation  ofr^. 

THEOREM  3.1:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that 
P{D‘^(0,y9)}  =  1.  At  iteration  h  =  1,2,. . .,  if  P{t^  <  h}  >  0  for  some  f  co*  e  G,  then 

E[x^|T^<h]  =  h-i;  [l-fl  [l-r(j,)ff)]]/[l-n  [i-rQM] 

T=i  y=i  7=1 

Proof:  See  Orosz  and  Jacobson  (2002a). 

The  expression  in  Theorem  3.1  clearly  shows  that  E[x^  |  x^  <  h]  is  no  greater  than  h.  Theorem  3.2 
provides  an  expression  for  the  conditional  variance  of  x^. 

Theorem  3.2:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that  P{D‘^(0,y5)} 
=  1.  At  iteration  h  =  1,2,.. .,  if  P{xy3  <  h}  >  0  for  some  f>f{co*),  co*  e  G,  then 

Var[x^|x^<h]  =  (h+1)^-E  [(2x+l)(l-n  [1 -r(j,;ff)])  /  [1- ft  [l-rG,i^]]] 

r=0  7=1  7=1 

-[h-E[l-fl  [l-r(j/)]]/[l-ft  [l-r(j,m^ 

T=I  7=1  7=1 

Proof:  See  Orosz  and  Jacobson  (2002a). 

Theorem  3.3  provides  expressions  for  the  conditional  expectation  of  x^  given  that  a  larger  or  smaller 
p  value  has  been  reached. 

Theorem  3.3:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that  P{D‘^(0,yS2)} 
=  1.  At  iteration  h  =  1,2,...,  assume  that  P{xy?  <  h}>0  for  some  >  foo*)  and  >f{a>*),  co*  e  G,  with 
Pi  <  P2.  Then  at  iteration  h, 

(a)  E[x^,|x^<h]=i  [Tr(x,yff,)n  [l-ra,A)]  /  [1 -ft  [l-r0,y52)]]] 

T=1  7=1  M 

+00  h 

+  i;  [x  P{x^,  =  X,  x^  <  h}  /  [1  -n  [l-r(j,y52)]]] 

r=/i+l  j=\ 

(b)  E[x^|x^,<h]  =  (h+l)-i:  [P{x^,<x,X7,i<h}/[l-ft  [l-r(j,M] 

r=0  7=1 

Proof:  See  Orosz  and  Jacobson  (2002a). 

The  expression  for  E[x^i  |  x^2  ^  ^,P\  <  1^2  in  Theorem  3.3  provide  measures  for  assessing  the  finite¬ 
time  performance  of  a  GHC  algorithm.  Note  that  one  difficulty  with  these  expressions  is  that  they 
contain  infinite  summations,  hence  are  impractical  to  compute. 

Theorem  3.4  provides  expressions  for  the  conditional  variance  of  xp  given  that  a  larger  or  smaller  p 
value  has  been  reached. 

Theorem  3.4:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that  P{D‘^(0,y?2)} 
=  1 .  At  iteration  h  =  1,2,. . .,  assume  that  P  {Xy3<h}  >  0  for  some  >f(jo*)  and  P2  ^  Gr,  with  Pi 

<  Pi.  Then  at  iteration  h, 

(a)  Var[x^,|x^2<h]=i:  [x^r(x,A)  ft  [l-r(j  A)]  /  [1  "ft  [l-r(j,A)]]] 

i-=i  7=1  7=1 

+  E  [t^P{T^i=T,T^<h}/[l-ft  [l-rO,A)]]] 

r=/t+l  j=\ 
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-[i  [^r(T,y9i)n  [i-raySi)]/[i-n  [i-ra/2)]]] 

r=l  M  M 

+  E  [TP{T^,=T,T^<h}/[l-n  [l-r(j,M]]f 

T=h+\  j=l 

(b)  Var[T^|T^i<h]  =  (h+l)'-i:  [(2T+l)P{T^<T,T^,<h}/[l-n  [l-r(jA)]]] 

T=0  j=\ 

-[(h+i)-i  [P{T^<T,T^i<h}/[i-n  [i-r(j,A)]]]f 

r=0  >1 

Proof:  See  Orosz  and  Jacobson  (2002a). 

These  results  describe  how  the  one-step  yff-acceptable  probability  and  the  random  variable  can  be 
used  to  obtain  finite-time  performance  measures  for  GHC  algorithms.  The  main  difficulty  with  these 
measures  is  that  they  are  difficult  to  compute  or  estimate.  The  following  description  and  results  show 
how  these  measures  can  be  computed  for  a  specific  GHC  algorithm,  the  S^A  algorithm.  Similar  results 
for  cyclical  simulation  annealing  (CSA)  are  reported  in  Orosz  and  Jacobson  (2002b). 

To  obtain  a  run  length  analysis  of  S^A,  properties  of  the  S^A  algorithms  can  be  exploited  to  obtain 
upper  and  lower  bounds  for  E[Tyj].  To  obtain  these  bounds  for  a  GHC  algorithm,  the  total  number  of 
iterations  executed  can  be  decomposed  into  a  set  of  cycles,  each  of  length  h,  where  the  cycle  number  of 
iteration  t  is  denoted  by 

jh(x)  =  rT/h1,  (6) 

the  ceiling  function  for  x  /  h  (i.e.,  the  smallest  integer  greater  than  or  equal  to  x  /  h).  To  obtain  the  cycle 
in  which  the  algorithm  visits  any  element  in  the  set  of  ^-acceptable  solutions  for  the  first  time,  define  the 
random  variable 

Jh(xy?)  =  min{jh(T)  >  1  :  jh(T)h  >y,T=  1,2,...}  =  fxp  /  h1 .  (7) 

Figure  3.1  depicts  the  relationship  between  jh(T)  and  t,  where  jh(x)  =  j  for  all  values  of  x  between  (j-l)h+l 
and  jh.  Therefore,  if  the  set  of  ;5-acceptable  solutions  is  visited  for  the  first  time  between  iterations  (jh(x)- 
l)h+l  and  jh(x)h,  then  Jh(T^)  =  jh(x). 


XfrX 


0  h  2h  Gh(x)-l)h  jh(x)h 
Figure  3.1:  Relationship  between  jh(x)  and  xp 

Lemma  3.3  summarizes  relationships  that  exist  between  the  events  {x^=  x}  and  {Jh(x^)  =  jh(x)}. 

Lemma  3.3:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that  P{D°(0,jff)}  = 
1 .  Suppose  that  the  algorithm  is  executed  for  x  >  0  iterations.  Then  for  cycle  length  0  <  h  <  x, 

(a)  P{x^=x}<P{Jh(x^)=jh(x)}, 

(b)  P{Jh(x^)  >  jh(x)}  =  P{x^  >  jh(x)h}, 

(c)  P{Jh(x^)  >  jh(x)}  <  P{x^  >  jh(x)h  -  j}  for  all  j  >  0, 

(d)  P  { Jh(x^)  >  jh(x)}  ^  P  {xys  >  jh(x)h  +  j  }  for  all  j  >  0, 

(e)  P{Jh(x^)  =  jh(x)}  =  P{x^  <  jh(x)h}  -  P{x^  <  (jh(x)-l)h}, 

(f)  P{x^  >  jh(x)h}  <  P{x^>  x}  <  P{x^  >  (jh(x)-l)h}, 

Proof:  See  Orosz  and  Jacobson  (2002a). 

Lemma  3 .4  obtains  an  expression  for  the  expectation  of  Jh(x^)  for  GHC  algorithms  in  terms  of  the  one- 
step  yff-acceptable  probability. 
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Lemma  3.4:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that  P{D‘'(0,)8)}  = 
1.  Then 

+00  jh 

E[JhM=  E  [n  [l-r(i,yff)]]. 

7=0  !=1 

Proof:  See  Orosz  and  Jacobson  (2002a). 

Theorem  3.5  uses  the  random  variable  Jh(t^)  to  obtain  upper  and  lower  bounds  for  E[Tyj],  the  expected 
number  of  iterations  to  visit  the  set  of  jff-acceptable  solutions  for  the  first  time. 

Theorem  3.5:  Consider  a  GHC  algorithm  execution  with  initial  solution  generated  such  that  P{D‘'(0,y5)} 
=  1.  Then 

1  +  h  [E[h(y)]  -1]  <  E[t;,]  <  1  +  h  [E[Jh(T^)]]  (8) 

Proof:  See  Orosz  and  Jacobson  (2002a). 

Theorem  3.5  provides  upper  and  lower  bounds  for  E[Ty9]  that  are  functions  of  h  and  [E[Jh(T^)].  The 
problem  in  applying  these  bounds  is  that  they  contain  infinite  summations,  hence  are  not  easily 
computable.  This  problem  is  difficult  to  overcome  for  GHC  algorithms  in  general.  However,  the 
cyclical  nature  of  S^A  can  be  exploited  to  circumvent  this  problem.  In  particular,  if  the  random  variable 
Jh('t^)  can  be  modeled  as  a  geometric  distribution  with  parameter  P{t^  <  h},  then  E[Jh(T^)]  =  l/Pfi:^  ^  h}, 
and  the  bounds  in  Theorem  3.5  simplify  to 

1  +  h  [P{t^ >  h}  /  P{t^  <  h}]  <  E[t^]  <  1  +  h  [1  /  P{x^  <  h}]. 

From  Lemma  3.2,  P{r^  <  h}  can  be  estimated  using  a  finite  length  execution.  Therefore,  E[Jh(T^)]  for  S^A 
can  be  estimated  using  information  obtained  from  replicating  each  algorithm’s  performance  over  a  single 
cycle  (i.e.,  over  the  first  h  iterations).  To  see  this,  for  cycle  m,  the  probability  that  the  algorithm  visits  the 
set  of  y9-acceptable  solutions  for  the  first  time  in  cycle  m  is 

=  m}  =  P{(m-l)h  +1  <  <  mh}.  (9) 

Lemma  3.5  uses  this  expression  to  show  that  if  the  one-step  y9-acceptable  probability  is  constant  over  all 
iterations  (i.e.,  r{\,P)  =  r  for  all  i=l,2,...),  then  Jh(xys)  is  a  geometric  random  variable  with  parameter  P{x;j 
<h}  =  l-(l-r)^ 

Lemma  3.5:  Consider  a  S^A  algorithm  execution  with  initial  solution  generated  such  that  P  {0*^(0 ,y9)}  =  1 
and  cycle  length  h  >  0.  If  T{i,P)  =  r  for  i  =  1,2,  ...,  then  Jh(T/s)  is  a  geometric  random  variable  with 
parameter  P{x^  <  h}=  1-  (l-r)’’. 

Proof:  See  Orosz  and  Jacobson  (2002a). 

For  S^A  algorithms,  since  the  cooling  schedule  is  fixed  at  all  iterations,  then  the  probability  of 
accepting  up  hill  moves  is  independent  of  the  iteration  number.  Therefore,  though  it  is  likely  to  be  very 
difficult  to  formally  prove  this  result,  it  may  be  reasonable  to  assume  that  the  values  are  constant 
for  all  i  =  1,2,...,  hence  the  result  in  Lemma  3.5  applies.  Theorem  3.6  provides  an  upper  and  lower 
bound  on  the  expected  number  of  iterations  to  reach  the  set  of  y9-acceptable  solutions,  under  this 
assumption. 

Theorem  3.6:  Consider  a  S^A  algorithm  execution  with  initial  solution  generated  such  that  P{D°(0,yff)}  = 
1 .  If  Jh(x^)  is  a  geometric  random  variable  with  parameter  P  (x^  <  h} ,  then 

1  +  h  [P{x^  >  h}/  P{x^<  h}]  <  E[x^]  <  1  +  h  / P{x^<  h}. 

Proof:  See  Orosz  and  Jacobson  (2002a). 

Theorem  3.7:  Consider  a  S^A  algorithm  execution  with  initial  solution  generated  such  that 
P {D‘’(0,y5)}=l .  If  Jh(x^)  is  a  geometric  random  variable  with  parameter  P (x^  <  h} ,  then 

E[xp  I  xp <  h]  P{xp <  h}  +  [1  +  h /  P{xp <  h}]  P{xp >  h}  <  E[x^] 
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<  E[t:p  I  Tp <  h]  P{Tp<  h}  +  [h  +  1  +  h  /  P{Tp <  h}]  P{Tp>  h} 
and 

Var[T^]  <  h'(  [P{t;,>  h}+l]  /  P{x^<  h}')  -  (E[Tp  |  xp <  h]  P{Tp ^  h}  +  [1  +  h  / P{Tp <  h}]  P{xp >  h}f. 
Proof:  See  Orosz  and  Jacobson  (2002a). 

Theorem  3.8  shows  that  the  upper  and  lower  bounds  for  E[xp]  presented  in  Theorem  3.7  are  tighter 
than  the  upper  and  lower  bounds  for  E[xp]  presented  in  Theorem  3.6. 

Theorem  3.8:  Consider  a  S^A  algorithm  execution  with  initial  solution  generated  such  that  P  {0*^(0 ,y?)}  = 
1.  If  Jh(x^)  is  a  geometric  random  variable  with  parameter  P{x^  <  h},  then 

E[xp  I  xp<  h]  P{xp <  h}  +  [h  +  1  +  h / P{xp^  h}]  P{xp>  h}  <  1  +  h /  P{x^<  h} 
and 

E[xp  I  xp<  h]  P{xp<  h}  +  [1  +  h  /  P{xp<  h}]  P{xp>  h}  >  1  +  h  [P{xp>  h}/ P{x^<  h}]. 

Proof:  See  Orosz  and  Jacobson  (2002a). 

Theorem  3.8  shows  that  the  upper  and  lower  bounds  for  E[xp]  in  Theorem  3.7  are  tighter  than  the 
bounds  given  in  Theorem  3.6.  Therefore,  under  the  assumption  that  Jh(Xys)  is  a  geometric  random  variable 
with  parameter  P{x^  <  h},  the  bounds  in  Theorem  3.7  should  be  used.  Note  that  under  this  assumption, 
these  results  can  only  be  applied  to  certain  GHC  algorithms,  such  as  S^A.  The  following  description 
exploits  these  results  for  S^A  algorithms  to  show  how  the  bounds  in  Theorem  3.7  can  be  computed. 

Computational  results  with  the  traveling  salesman  problem  to  illustrate  how  the  upper  and  lower 
bounds  for  E[x^]  in  Theorem  3.7  can  be  computed.  The  traveling  salesman  (optimization)  problem  (TSP) 
is  a  well-studied  NP-hard  discrete  optimization  problem  (Lawler  et  al.  1985).  The  diversity  of 
applications  for  the  TSP  makes  it  a  frequent  choice  for  testing  and  evaluating  the  efficiency  and 
effectiveness  of  algorithms  and  heuristics  for  intractable  discrete  optimization  problems.  Traditional 
applications  of  the  TSP  can  be  found  in  numerous  domains,  including  vehicle  routing  and  scheduling 
problems.  More  recently,  it  has  been  applied  to  the  manipulation  of  robotics  (Balaguer  et  al.  2000),  the 
cutting  of  industrial  components  (Foerster  and  Wascher  1998),  and  circuit  board  design  (Kobayashi  et  al. 
1999).  A  search  and  rescue  military  application  can  also  be  modeled  with  the  traveling  salesman 
problem  (Henderson,  Vaughan,  and  Jacobson  2003). 

To  apply  a  local  search  algorithm  to  an  instance  of  the  TSP,  a  neighborhood  function  must  be 
defined.  There  are  numerous  neighborhood  functions  that  have  been  devised  for  the  TSP.  One  such 
neighborhood  function  is  the  4-change  method.  The  4-change  method  moves  between  solutions  by  the 
exchange  of  four  edges.  There  are  several  methods  that  accomplish  the  4-change  neighborhood  function. 
One  such  method  is  the  city-exchange  method.  The  city  exchange  method  defines  a  neighbor  of  a  given 
cycle  by  exchanging  two  cities  in  the  cycle.  For  example,  for  a  set  of  seven  cities  C  =  {ci,  Ca,  C3,  C4,  Cs, 
Ce,  C7}  with  solution  space  Q.  =  {coi,  coa,  ...,  ®36o},  if  the  current  solution  is  00  =  (c’l,  c’a,  c’3,  c’4,  c’5,  c’e, 
c’7),  then  by  exchanging  cities  c’l  and  c’5,  a  city  exchange  neighboring  solution  is  co’  =  (c’5,  c’a,  c’3,  c’4, 
c’i,c’6,  c’  7),  with  four  edges  exchanged. 

Another  commonly  used  neighborhood  function  for  the  TSP  is  the  2-Opt  neighborhood  function 
(Croes  1958).  The  2-Opt  neighborhood  function  moves  from  one  solution  to  another  solution  by  the 
exchange  of  two  edges.  For  example,  consider  the  finite  set  of  seven  cities  C  =  {ci,  Ca,  C3,  C4,  C5,  Ce,  C7}, 
and  the  corresponding  solution  space  Q  =  {®i,  coa,  ...,  (Oaeo).  If  the  current  solution  is  ro  =  (c’l,  c’a,  c’3, 
c’4,  c’5,  c’s,  c’7),  then  one  possible  neighboring  solution  is  co’  =  (c’l,  c’5,  c’4,  c’3,  c’a,  c’e,  c’7),  which  is 
obtained  by  reversing  the  sequence  of  cities  (c’a,  c’3,  c’4,  c’5).  The  2-Opt  neighborhood  function  is  a 
specific  version  of  a  more  general  neighborhood  function  termed  X.-Opt  (Helsgaun  2000),  where  in  each 
move  from  one  solution  to  another  solution,  X  edges  are  exchanged.  The  X,-opt  neighborhood  function  is 
based  on  the  concept  that  a  cycle  is  considered  ^.-optimal  if  it  is  impossible  to  obtain  a  lower  objective 
function  value  by  the  exchange  of  any  X  edges.  From  this  definition,  as  X  increases,  the  resulting  local 
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search  algorithm  is  more  likely  to  find  an  optimal  solution.  Unfortunately,  the  number  of  operations  to 
test  all  ^.-exchanges  also  increases  as  the  number  of  cities  increases.  Therefore,  local  search  algorithms 
that  use  the  X.-C)pt  neighborhood  function  typically  restrict  the  value  for  X  to  be  two  or  three  (Helsgaun 
2000). 

The  Lin-Kemighan  neighborhood  function  (Lin  and  Kergnighan  1973)  operates  by  determining  the  X 
value  in  the  X-opt  neighborhood  function  that  achieves  the  best  compromise  between  algorithm  run-time 
and  quality  of  solution.  At  each  iteration  of  the  algorithm,  a  series  of  tests  is  performed  to  determine  the 
optimal  number  of  X.-exchanges  that  should  be  considered.  This  method  is  commonly  considered  the 
most  effective  neighborhood  function  as  well  as  the  most  difficult  to  design.  Note  that  since  the  purpose 
of  the  computational  experiments  is  not  to  determine  an  optimal  neighborhood  function  for  local  search 
algorithms  applied  to  the  TSP,  but  rather,  to  illustrate  the  upper  and  lower  bounds  derived  in  Section  4, 
only  the  city  exchange  neighborhood  function  is  used. 

The  S^A  algorithm  was  applied  to  a  randomly  generated  one  hundred  city  instance  of  the  TSP.  The 
distance  matrix  D(Ci,Cj)  between  each  pair  of  cities  Ci,Cj  e  C  was  generated  by  randomly  placing  the  cities 
on  a  10  X  10  Cartesian  grid  and  using  a  two-dimensional  array  to  store  the  (X,Y)  coordinates  of  each  city 
(each  coordinate  was  independently  generated  uniform  [0,10]).  The  distance  between  cities  was  then 
computed  and  placed  in  the  distance  matrix.  Note  that  by  design,  the  distance  matrix  is  symmetric  (i.e., 
D(Ci,Cj)  =  D(Cj,Ci)  for  each  pair  of  cities  Ci,Cj  e  C). 

The  S^A  algorithm  was  applied  with  six  different  (fixed)  temperatures  (T  =  0,  1,5,  20,  50,  100)  and 
four  different  cycle  lengths  (h  =  500,  2000,  5000,  and  2,000,000),  resulting  in  a  total  of  twenty-four 
different  parameter  settings  (hence  algorithm  executions).  At  each  iteration  of  the  S^A  algorithm,  a 
neighbor  of  the  current  cycle  was  generated  using  the  city  exchange  neighborhood  function,  where  the 
two  cities  selected  for  exchange  were  uniformly  generated  across  each  of  the  cities  (note  that 
experiments  using  the  2-opt  neighborhood  function  were  also  performed,  but  the  results  obtained  were 
inferior,  as  measured  by  the  values  for  E[Tp],  to  those  obtained  using  the  city  exchange  neighborhood). 
For  each  S^A  algorithm  execution,  an  initial  cycle  was  randomly  selected  among  all  possible  Hamiltonian 
cycles.  Based  on  all  the  S^A  algorithm  executions  and  replications,  the  smallest  Hamiltonian  cycle 
obtained  had  total  distance  99.4,  which  serves  as  an  upper  bound  on  the  optimal  Hamiltonian  cycle 
distance.  One  thousand  replications  were  executed  for  each  of  the  parameter  settings  to  estimate  the  one- 
step  ;5-acceptable  probability  at  each  iteration.  For  each  of  these  replications,  the  same  problem  instance 
and  parameter  settings  were  used,  with  a  different  randomly  generated  initial  cycle.  The  average  initial 
solution  over  the  24,000  runs  performed  (twenty-four  parameter  settings  with  one  thousand  runs  each) 
had  a  distance  of  approximately  479.  Each  of  the  twenty-four  experiments  was  executed  on  a  450  MHz 
Dell  PC.  For  the  cycle  lengths  h  =  500,  2000,  5000,  and  2,000,000  the  PC  completed  all  six  of  the  fixed 
temperature  settings  in  approximately  24  seconds,  1.6  minutes,  4  minutes,  and  25  hours,  respectively. 
Different  y9-acceptable  values  are  reported  for  each  parameter  setting  of  T  and  h.  Note  that  the  bounds 
for  only  y^-acceptable  values  that  are  reached  over  all  1000  replications  for  the  h  =  2,000,000  executions 
are  reported;  this  allows  the  results  for  the  lower  and  upper  bounds  from  Theorem  3.7  for  smaller  values 
of  h  to  be  validated. 

Orosz  and  Jacobson  (2002a,b)  report  extensive  computational  results  with  the  lower  and  upper 
bounds  for  E[Tp]  in  Theorem  3.7  for  the  S^A  runs  with  cycle  length  h  =  500,  2000,  for  fixed  temperature 
values  T  =  0,  1,5,  20,  50,  and  100.  These  results  suggest  that  the  lower  and  upper  bound  estimators  are 
most  effective  for  larger  temperatures.  Therefore,  as  the  temperature  parameter  increases,  the  accuracy 
for  predicting  a  lower  and  upper  bound  on  E[Tp]  improves.  This  is  not  surprising  since  as  the  temperature 
parameter  increases,  the  S^A  algorithm  traverses  the  solution  space  more  randomly,  hence  Jh(x^)  is  more 
likely  to  be  modeled  as  a  geometric  random  variable  with  parameter  P{t^  <  h}.  Furthermore,  since  the 
lower  and  upper  bounds  and  the  values  for  E[Tp]  are  estimates,  then  the  upper  bound  estimates  may 
actually  be  less  than  the  estimates  for  E[Tp].  Conversely,  the  lower  bound  estimates  may  actually  be 
greater  than  the  estimates  for  E[Tp].  This  error  in  the  lower  and  upper  bound  estimates  appears  to  occur 
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when  P{Tp  <  h}  <  0.20.  This  may  also  be  due  to  the  geometric  distribution  assumption  on  IhCi:;!)  not  being 
valid. 

The  computational  results  reported  in  Orosz  and  Jacobson  (2002a, b)  with  a  one  hundred  city  instance 
of  the  TSP  validate  the  lower  and  upper  bound  expressions  for  E[Tp]  in  Theorem  3.7.  These  results 
suggest  that  the  lower  and  upper  bound  expressions  for  E[Tp]  in  Theorem  3.7  provide  reasonable 
measures  at  higher  constant  temperature  parameter  settings.  In  particular,  if  similar  values  of  P{Tp  <  h} 
are  compared  across  different  constant  temperature  parameter  settings,  T,  for  the  same  cycle  length,  h, 
then  the  upper  and  lower  bound  estimates  for  E[Tp]  are  most  accurate  for  higher  values  of  T.  Moreover, 
the  estimate  for  E[Tp]  falls  between  the  upper  and  lower  bound  estimates  for  E[xp]  for  approximately  40% 
of  the  data  reported  over  all  six  tables  in  the  upper  data  set  where  P{Tp  <  h}  <  1.  However,  of  the 
instances  where  the  estimate  for  E['rp]  does  not  fall  between  the  upper  and  lower  bound  estimates  for 
E[Tp],  the  estimates  at  higher  constant  temperature  parameter  settings  are  much  closer  to  the  upper  and 
lower  bound  estimate  range  than  the  estimates  at  lower  constant  temperature  parameter  settings.  These 
computational  results  also  suggest  that  the  upper  and  lower  bound  estimates  for  E[Tp]  have  a  small 
standardized  error  at  higher  constant  temperature  parameter  settings,  hence  provide  useful  information  as 
estimates  for  E[Tp].  Lastly,  these  results  suggest  that  the  lower  and  upper  boimd  expressions  for  E[xp]  in 
Theorem  3.7  provide  better  measures  when  estimated  with  longer  cycle  lengths.  This  relationship  is 
reasonable,  since  as  the  cycle  length  increases,  the  value  for  P{xp<h}  for  each  value  of  P  also  increases. 

4.  Construction  Site  Leveling  Problem 

An  interesting  application  of  the  research  results  obtain  has  been  a  construction  site  leveling  problem,  the 
details  of  which  are  reported  in  Henderson  et  al.  (2003).  To  describe  the  results  obtained,  the  following 
background  information  is  provided.  Heavy  engineering  and  construction  projects  often  require  terrain 
modifications  that  involve  moving  large  amounts  of  earth  from  one  area  to  another  area  of  a  construction 
site.  Excavating  earth  from  cut  (surplus)  locations  and  hauling  it  for  deposit  into  fill  (deficit)  locations 
requires  earthmoving  vehicles  that  are  expensive  to  both  operate  and  maintain,  hence  often  constitute  a 
large  portion  of  a  construction  project’s  budget  (see  Bartholomew  2000).  Ideally,  planners  must  develop 
a  strategy  for  contouring  terrain  that  minimizes  the  total  distance  traveled  by  earthmoving  vehicles 
between  cut  and  fill  locations.  Minimizing  the  total  distance  traveled  results  in  an  overall  savings  to  the 
construction  project. 

The  Shortest  Route  Cut  and  Fill  Problem  (SRCFP)  seeks  to  find  a  route  (beginning  and  ending  at  the 
same  cut  location)  for  a  single  earthmoving  vehicle  that  minimizes  the  total  distance  traveled  between  cut 
and  fill  locations.  Complete  enumeration  of  all  possible  solutions  to  the  SRCFP  would  take  a  prohibitive 
amount  of  time.  Therefore  it  is  necessary  to  construct  efficient  and  effective  optimization  algorithms  to 
identify  optimal/near-optimal  routes  for  construction  project  planners. 

One  approach  to  finding  solutions  to  the  SRCFP  is  to  apply  a  greedy  algorithm.  A  greedy  algorithm 
applies  a  myopic  strategy  that  begins  at  an  arbitrary  cut  location,  moves  to  the  nearest  fill  location,  and 
then  moves  from  that  fill  location  to  the  nearest  cut  location.  This  process  of  moving  to  the  nearest 
available  location  continues  until  the  construction  project  site  is  leveled  (i.e.,  all  the  fill  locations  have 
been  filled  with  earth  fi'om  the  cut  locations).  This  report  summarizes  results  obtained  using  this  greedy 
algorithm.  The  cost  associated  with  the  solution  found  using  the  greedy  algorithm  provides  an  upper 
bound  for  the  optimal  solution  for  the  problem.  Local  search  algorithms  have  also  be  used  to  find  near- 
optimal  solutions  to  the  SRCFP.  This  report  also  applies  simulated  annealing  (one  such  local  search 
algorithm)  to  the  SRCFP.  The  report  is  organized  as  follows:  First,  definitions  that  are  needed  to 
describe  the  SRCFP  and  a  brief  description  of  the  simulated  annealing  algorithm  are  presented.  The 
SRCFP  is  then  formally  stated  as  an  NP-Hard  discrete  optimization  problem.  Computational  results  for 
ninety  instances  of  the  SRCFP  with  simulated  annealing  algorithms  are  reported.  Upper  and  lower 
bounds  for  the  optimal  objective  function  value  for  the  SRCFP  instances  are  also  given. 

To  describe  the  SRCFP,  several  definitions  and  terminology  are  needed.  Define  a  construction 
project  site  as  a  plot  of  land  with  existing  contours  that  can  be  modified  to  facilitate  construction  of  an 
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object  (e.g.,  road,  building,  runway).  A  field  survey  of  the  construction  project  site  is  typically 
preformed  before  any  work  begins  on  the  construction  project.  This  field  survey  results  in  a  uniform  grid 
that  defines  the  construction  project  site  by  coordinate  pairs.  Define  locations  on  the  construction 
project  site  as  coordinate  pairs  resulting  from  a  field  survey  of  the  construction  project  site.  A  cut 
location  is  a  location  on  the  construction  project  site  with  excess  (surplus)  earth  and  a  fill  location  is  a 
location  on  the  construction  project  site  that  requires  earth  (deficit).  The  final  grade  is  achieved  when  all 
excess  earth  at  cut  locations  is  removed  and  all  earth  deficits  at  fill  locations  are  filled.  Note  that  not  all 
locations  have  surplus  or  deficit  earth  (i.e.,  a  particular  location  elevation  may  be  equal  to  the  final 
grade).  Define  a  unit  load  as  the  volume  of  earth  that  an  individual  vehicle  can  carry  in  one  trip  (e.g.,  a 
sixteen  cubic  yard  scraper  versus  a  thirty  cubic  yard  dump  truck).  Define  the  number  of  unit  cuts  as  the 
number  of  unit  loads  available  (excess)  at  a  cut  location  and  the  number  of  unit  fills  as  the  number  of  unit 
loads  required  (deficit)  at  a  fill  location.  The  construetion  project  site  contains  m  unit  cuts  and  n  unit 
fills.  To  guarantee  that  a  feasible  solution  exists  for  the  SRCFP,  assume  that  adequate  offsite  earth  is 
available  to  ensure  that  the  total  deficit  can  be  satisfied  and/or  that  a  disposal  site  is  available.  This 
assumption  implies  m  =  n,  since  any  excess  cut  (fill)  location  unit  loads  can  be  matched  with  dummy  fill 
(cut)  location  unit  loads.  A  route  is  the  path  (beginning  and  ending  at  the  same  cut  location)  that  an 
earthmoving  vehicle  follows  to  visit  every  unit  cut  and  every  unit  fill  location  exactly  once,  alternating 
between  the  cut  and  fill  locations,  until  the  terrain  modifications  are  complete  then  returning  to  initial 
location.  The  total  haul  distance  is  the  length  of  a  given  route. 

The  objective  of  the  SRCFP  is  to  find  a  route  that  minimizes  the  total  haul  distance  traveled  by  a 
vehicle  (i.e.,  a  Hamiltonian  circuit  with  alternating  unit  cut  and  unit  fill  locations  that  transforms  a  given 
piece  of  terrain  into  the  final  grade).  Since  the  SRCFP  seeks  a  route  that  minimizes  the  total  haul 
distance  while  visiting  every  location  exactly  once,  the  SRCFP  is  a  special  case  of  the  symmetric 
traveling  salesman  problem.  The  SRCFP  can  also  be  related  to  the  capacitated  traveling  salesman 
problem  with  pickups  and  deliveries  (CTSPPD;  see  Anily  and  Bramel  1999),  since  it  consists  of  two 
location  types  and  uses  a  vehicle  with  limited  capacity.  Therefore,  SRCFP  is  a  special  case  of  the 
swapping  problem  where  the  vehicle  capacity  is  equal  to  one  (Anily  and  Hassin  1992).  Chalasani  and 
Motwani  (1995)  presents  a  2-approximation  algorithm  to  address  a  special  case  of  the  swapping  problem 
with  two  product  types,  which  is  equivalent  to  a  special  case  of  the  SRCFP  with  vehicle  capacity  one. 

Note  that  the  SRCFP  can  be  formulated  as  a  transportation  problem  by  relaxing  the  requirement  that  a 
route  must  alternate  between  unit  cut  locations  and  unit  fill  locations,  and  by  modifying  the  objective 
function  to  minimize  the  cost  of  the  distance  traveled  between  locations.  The  transportation  problem 
optimally  allocates  resources  by  minimizing  the  distance  traveled  between  cut  and  fill  locations;  however, 
it  does  not  take  into  account  the  cost  of  traveling  from  each  fill  location  to  the  next  cut  location. 
Therefore,  relaxing  the  SRCFP  to  the  transportation  problem  results  in  the  concept  of  routes  being  lost. 

Local  search  algorithms  require  a  neighborhood  function  at  each  solution  in  the  solution  space. 
Unfortunately,  the  very  design  of  local  search  algorithms  means  that  they  often  reach  a  local  optimum 
(for  a  given  neighborhood  function),  and  may  not  be  able  to  escape  this  local  optimum  to  continue 
searching  for  global  optima  (e.g.,  deterministic  local  search;  see  Aarts  and  Lenstra  1997).  Simulated 
annealing  (Henderson,  Jacobson,  and  Johnson  2003)  is  a  local  search  algorithm  used  to  address  hard 
discrete  optimization  problems.  Simulated  annealing  allows  for  the  escape  from  local  optima,  with  the 
possibility  of  reaching  a  global  optimum,  by  allowing  uphill  moves.  To  describe  the  implementation  of 
simulated  annealing,  the  following  definitions  are  needed.  Let  Q  be  the  solution  space  (i.e.,  the  set  of  all 
possible  solutions).  Let  f:  Q  9?  be  the  objective  function  defined  on  the  solution  space.  The  goal  is  to 
find  a  global  minima,  ©*,  (i.e.,  oo*  e  Q  such  that  f(ft))  >  f(o)*)  for  all  ro  e  Q).  Note  that  the  objective 
function  is  assumed  to  be  bounded  to  ensure  that  co*  exists.  Define  ri(co)  to  be  the  neighborhood  function 
for  ©  e  Q.  Therefore,  associated  with  every  solution,  ©  e  Q,  are  neighboring  solutions,  ti(©),  that  can  be 
reached  in  a  single  iteration  of  the  simulated  annealing  algorithm.  Q.  Define  gij(k)  to  be  the  generation 
probability  function  for  the  neighborhood  function  q,  where  the  probability  that  ©jeq(©i)  is  generated 
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during  iteration  k,  is  gij(k).  Simulated  annealing  starts  with  an  initial  solution  oo  e  Q.  A  neighboring 
solution  co’  6  ri((o)  is  then  generated.  If  f(co’)  >  f(co),  then  co’  is  accepted  as  the  current  solution  with 
probability  e'^  ^ were  T(t)  is  a  temperature  parameter  that  is  typically  non-increasing  at  each 

iteration,  t. 

The  SRCFP  can  be  formulated  as  a  discrete  optimization  problem  with  the  objective  of  minimizing 
the  total  haul  distance  traveled  by  an  earthmoving  vehicle  between  cut  and  fill  locations.  The  number  of 
visits  to  a  location  depends  on  the  number  of  unit  loads  necessary  to  level  the  location.  Since  each  unit  of 
earth  represents  a  unit  load,  the  number  of  visits  to  each  location  is  the  amount  of  excess  (or  deficit)  earth 
at  that  location. 

To  model  the  SRCFP  as  a  discrete  optimization  problem,  several  definitions  are  needed.  Define  G  = 

{gi,  g2, . gn}  to  be  a  set  of  n  locations.  Define  V  =  {vi,  V2,....,  Vn}  to  be  a  set  of  volumes  for  each 

location  in  G,  (i.e.,  the  number  of  unit  cuts  or  unit  fills).  Define  L  =  {//,  h, 4}  to  be  the  set  of  single 
unit  cut  and  unit  fill  locations,  that  correspond  to  locations  in  G,  where  gi  e  G  appears  in  L  exactly  |vj| 
times.  Therefore,  each  element  in  L  represents  a  single  (either  positive  or  negative)  volume  of  cut  or  fill, 

hence  X|v/|  =  k.  Define  =  {t  1, 1^2, .....  tm}  to  be  the  set  of  unit  cut  locations  and  L'  -  {Fi,  42, ....  ./"m) 

/=1 

to  be  the  set  of  unit  fill  locations,  where  L=V''^  L',\l\  =  k,  and  I I  =  I Z’  I  =  4  /  2  =  m.  Recall  that  the 
number  of  unit  cut  locations  equals  the  number  of  unit  fill  locations,  and  that  a  unit  cut  and  a  unit  fill 
location  cannot  occur  at  the  same  location  in  G.  Define  Cj  e  L*'  to  be  a  single  unit  cut  location,  i  = 
1,2,... ,m,  and  fi  e  L'  to  be  a  single  unit  fill  location,  i  =  1,2,... ,m.  Lastly,  let  R  =  (ci,  fi,...,  c^,  /„) 
represent  a  route,  defined  as  a  Hamiltonian  circuit  over  L  with  the  added  constraint  that  the  route  must 
alternate  between  unit  cut  locations  and  unit  fill  locations.  The  SRCFP  discrete  optimization  problem  is 
now  formally  stated. 

Shortest  Route  Cut  aud  Fill  Problem  (SRCFP) 

INSTANCE:  Given  a  set  of  m  unit  cut  locations,  L*=  {t  1X2,  ■■■Xm],  and  a  set  of  m  unit  fill  locations,  L' 
=  1,1^2,  ■■■Jm},  and  a  symmetric  matrix  D  containing  the  distance  between  the  cut  and  fill  locations. 

QUESTION:  Find  a  route  R  =  (cj,  fi,...,  f„),  c;  &  L*,\  =  1,2,. ..,m,  fie  L  ,i  =  1,2,. ..,m,  such  that 

m-\ 

giR)  =  E  [Picxfi^  +  D(4,Ci+i)]  +  D(c^r«)  +  D(fi^,Ci)  is  minimized. 

/=1 

The  SRCFP  includes  the  symmetric  traveling  salesman  (optimization)  problem  (STSP)  (Aarts  and 
Lenstra  1997)  as  a  special  case.  Therefore,  the  SRCFP  is  NP-Hard  (Garey  and  Johnson  1979).  Define  Q 
to  be  the  solution  space  of  routes  with  alternating  unit  cut  and  unit  fill  locations  (i.e.,  Q  =  {{cj,  fi,...,  Cm, 
fn),  Ci  G  L^,f\  G  Z",  i  =  1,2,. .  .,m}).  The  cardinality  of  Q  depends  on  the  number  of  unit  cut  and  unit  fill 
locations  for  the  given  problem  instance.  Lemma  4.1  provides  a  closed  form  expression  for  this  number, 
and  shows  that  the  size  of  the  solution  space  grows  exponentially  as  the  number  of  unit  cut  locations  and 
units  fill  locations  increase. 

Lemma  4.1:  Given  an  instance  of  the  SRCFP  with  m  unit  cut  and  m  unit  fill  locations,  the  size  of  the 
solution  space  is  |Q|  =  {m-\)\  m\  /  2. 

Proof:  See  Henderson  et  al.  (2003). 

Computational  results  were  reported  with  simulated  annealing  applied  to  ninety  randomly  generated 
instances  of  the  SRCFP.  Each  instance  of  the  SRCFP  corresponds  to  a  construction  project  site 
consisting  of  a  square  grid  with  25,  49,  or  81  locations  (i.e.,  a  construction  site  survey  resulting  in  5  x  5, 
7x7,  and  9x9  locations).  Thirty  instances  with  each  such  grid  size  are  reported.  The  distance  matrices 
are  constructed  with  one  unit  of  distance  between  each  location  in  the  horizontal  direction  and  one  unit  of 
distance  between  each  location  in  the  vertical  direction.  The  Euclidean  distance  between  locations  i  and 
7,  7  =  1,  2,  . .  .m,  is  denoted  by  Ay-  For  each  instance,  the  volume  of  earth  (i.e.,  number  of  unit  cuts  or 
unit  fills)  at  every  location  is  generated  as  a  discrete  uniform  [-3,  3]  random  variable.  Moreover,  each 
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problem  instance  includes  an  offsite  location  one  unit  of  distance  vertically  from  location  In.  This 
location  guarantees  the  assumption  that  adeqimte  offsite  earth  is  available  to  ensure  that  the  total  deficit 
can  be  satisfied  and/or  that  a  disposal  site  is  also  available,  hence  a  balanced  construction  project  site  is 
generated  (i.e.,  m  =  ri).  Note  that  the  random  volume  of  earth  generated  at  each  location  determines  the 
actual  size  of  the  problem  instance  (i.e.,  the  number  of  unit  cut  and  unit  fill  locations).  For  example,  the 
largest  9x9  locations  SRCFP  problem  instance  can  in  theory  include  as  many  as  243  total  unit  cut  and 
unit  fill  locations,  though  on  average,  such  problem  instances  will  contain  approximately  121.5  unit  cut 
and  unit  fill  locations.  From  Table  4.3,  all  thirty  randomly  generated  instances  of  the  SRCFP  for  this  size 
contained  between  100  and  134  total  unit  cut  and  unit  fill  locations,  hence  provide  large  SRCFP  problem 
instances  to  assess  the  effectiveness  of  simulated  annealing. 

For  simulated  annealing,  the  cooling  schedule  {T(t)}  is  updated  by  multiplying  the  previous 
temperature  parameter  by  an  increment  multiplier,  P,  where  0  <  P  <  1  (i.e.,  T(t)  =  pT(t-l)).  The  initial 
temperature,  T(0)  =  .2  (C  +  F)M  where  C  =  number  of  unit  cut  locations,  F  =  number  of  unit  fill 
locations  and  M  =  Max {!)(/;, /y),  ij  =  1,  2,. . .,  m}  with  P  =  .98  for  the  ninety  randomly  generated  instances 
of  the  SRCFP.  The  neighborhood  function  t|  implemented  for  all  experiments  is  defined  as  follows:  For 
all  routes  R  =  (cy,  /y,  C2,  fi,  •••Cm,  fm),  the  neighbors  of  R,  r|(R)  ,  are  defined  by  selecting  any  c’,c"  e  {cj, 
Cl,  ■■■,  Cm},  c'  ^  c",  and  reversing  the  sequence  of  locations  (both  cut  and  fill)  between  them.  Therefore 

fi(^)  .  R  (Cj,  J'l,  Cl,  Jl,“.Ci.j,  fi.i,  Cj,  fj.jCj.i,  ...,  fi+l^  Ci+j,  fi,  Cj,  fj,  Cj+j,  fj+i,  ...,  Cm,  fm}, 

for  some  i,j  =1,2,. . .,m,  i  <  j) . 

This  neighborhood  function  is  similar  to  2-opt  for  the  traveling  salesman  problem  (Aarts  and  Lenstra 
1997).  Note  that  since  the  distance  matrix  is  symmetric,  this  neighborhood  function  includes  all  routes 
where  two  fill  locations  are  selected  and  the  sequence  of  locations  between  them  are  reversed.  To 
generate  a  neighbor  at  each  iteration,  d  is  generated  uniformly  over  the  set  {ri,  fi,...,  Tn,}  and  c"  is 
generated  uniformly  over  the  set  {  (fi,  ^"2, . . .,  Tm}  \  (c'} } . 

Computational  results  with  simulated  annealing  are  reported  for  the  ninety  randomly  generated 
instances  of  the  SRCFP.  Simulated  aimealing  was  executed  with  x  =  500  and  N(t)  =  2{C+F)  for  t  =  1,  2, 
...,  X.  Fifty  replications  applying  simulated  annealing  to  each  problem  instance  were  executed.  The 
initial  solution  for  each  replication  was  randomly  generated  by  considering  two  sets:  one  of  unit  cut 
locations  and  one  of  unit  fill  locations.  These  two  sets  are  randomly  permuted  to  generate  the  resulting 

unit  cut  and  unit  fill  locations  for  the  initial  route.  The  initial  route,  R ,  was  then  reconstructed  by 
merging  the  permuted  unit  cut  and  unit  fill  sets  alternating  unit  cut  and  unit  fill  locations.  Initial 
solutions  for  each  replication  were  obtained  in  the  same  manner  by  randomly  permuting  the  first 
replication’s  initial  solution. 

The  mean  (p),  the  standard  deviation  (a),  and  the  minimum  and  maximum  of  the  objective  function 
values  over  the  fifty  replications  are  reported.  Tables  4. 1-4.3  report  results  for  the  fifty  replications  of 
the  simulated  annealing  algorithm  for  each  of  the  ninety  randomly  generated  instances  of  the  SRCFP. 
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Table  4.1 

Simulated  Annealing  Results 
(5  Unit  Cut  Locations  and  5  Unit  Fill  Locations) 

Problem 

iC+F) 

P 

a 

Minimum 

Maximum 

5x5  1 

36 

41.36 

0.48 

41.14 

42.37 

5x5  2 

38 

47.74 

0.32 

47.37 

48.02 

5x5  3 

36 

61.04 

0.08 

60.97 

61.40 

5x5  4 

38 

49.16 

0.12 

49.00 

49.25 

5x5  5 

36 

59.57 

0.62 

58,80 

60.73 

5x5  6 

46 

68.36 

0.59 

68.02 

70.13 

5x5  7 

30 

45,14 

0.60 

44.50 

5x5  8 

38 

64.48 

0.11 

64.47 

5x5  9 

40 

77.24 

0.00 

77.24 

5x5  10 

44 

49.08 

0.34 

48.55 

5x5  11 

38 

54.72 

mm 

54.65 

55.24 

5x5  12 

32 

43.71 

mm 

43.66 

44.25 

5x5  13 

40 

51.50 

0.68 

51.02 

5x5  14 

36 

60.88 

0.22 

60.78 

5x5  15 

40 

62.56 

0.32 

62.39 

63.19 

5x5  16 

44 

54.14 

0.46 

53.90 

55.31 

5x5  17 

34 

54.27 

0.30 

54,14 

54.90 

5x5  18 

34 

46.38 

0.23 

46.34 

47.51 

5x5  19 

24 

38.54 

0.46 

38.33 

39.51 

5x5  20 

26 

40.60 

40.60 

40.60 

5x5  21 

30 

37.14 

0.41 

36,96 

38.31 

5x5  22 

34 

50.26 

49.98 

50.68 

5x5  23 

44 

58.66 

0.40 

58.19 

59.74 

5x5  24 

36 

58.44 

58.09 

60.83 

5x5  25 

42 

66.86 

66,84 

68.06 

5x5  26 

42 

88.04 

0.04 

88.01 

88.21 

5x5  27 

38 

65.42 

0.05 

65.36 

65.47 

5x5  28 

32 

41.48 

0.32 

41.38 

42.44 

5x5  29 

42 

57.32 

0.19 

57.16 

57.78 

5x5  30 

48 

55.48 

0.47 

54.61 

55.86 

72 

78 

78 

Table  4.2 

Simulated  Annealing  Results 
7  Unit  Cut  Locations  and  7  Unit  Fill  Locations 


Problem  I  {C+F) _ |x _ ct  Minimum 

60  101.45  0.00  101.45 


I  106.34  0.89  105.11 

07.17  0.57  106.31 

30.00  0,07  129.90 

66  I  99.45  0.29  98.97 

05.32  0.87  104.42 

60.77  0.59  160.00 

63.89  0.39  163.21 

68  I  112.13  0.76  111.62 

80 _ 118.76  0,35  118.01 

80 _ 134.47  0.48  133.86 

76 _ 104.70  0.50  103.95 

62  85.50  1.20  84.07 


7x7  14 _ 70 _ 106.22  1.02  105.76 

7x7  15 _ 72 _ 129.37  0.72  128.403 

7x7  16 _ 64 _ 94.22  0.60  93.96 

_ 66  ^104.67  I  0.32 

7x7  18 _ 66 _ 124.96 

7x719  I  80  ~  102.39 

56 _ 115.50 

68 _ 98.23  I  0.67  _ 

68  95.40  0.70 


68 

108.68 

66 

107.42 

72 

95.95 

84 

140.17 

86 

117.90 

76 

115.20 

Maximum 

101.45 


08.01 
09.27 
30.13 
100.30 


07 

62 


164.65 


114.17 

119.59 

135.82 

106.19 

90.52 


110.82 


131.42 

96.33 


118.58 


00.32 

98.44 


0.52 

108.25 

111.54 

0.50 

106.83 

108.83 

0.71 

94.86 

97.04 

0.86 

139.54 

141.79 

.62 

117.18 

119.21 

,05 

113.90 

119.60 

ESI 

03 

99.24 

123.57 

Table  4.3 

Simulated  Annealing  Results 
(9  Unit  Cut  Locations  and  9  Unit  Fill  Locations) 


Problem 

(C+F) 

p 

a 

Minimum 

Maximum 

9x9  1 

122 

245.56 

0.68 

244.62 

246.87 

9x9  2 

132 

173.03 

1.73 

170.17 

9x9  3 

124 

203.55 

0.70 

202.25 

9x9  4 

116 

178.91 

1.62 

177.00 

185.75 

9x9  5 

118 

215.41 

214.34 

217.41 

9x9  6 

132 

212.51 

1.12 

210.53 

215.62 

9x9  7 

126 

204.15 

0.98 

202.53 

206.83 

9x9  8 

134 

mmm 

186.83 

190.84 

9x9  9 

116 

mssm 

191.87 

195.49 

9x9  10 

224.54 

0.57 

223.66 

225.35 

9x9  11 

219.78 

0.99 

217.93 

223.31 

9x9  12 

116 

243.29 

0.47 

242.49 

244.34 

9x9  13 

116 

172.06 

0.96 

170.40 

174.16 

9x9  14 

112 

169.44 

1.03 

168.10 

172.85 

9x9  15 

122 

217.55 

0.84 

216.33 

219.90 

9x9  16 

112 

156.03 

MKM 

153.78 

160.17 

9x9  17 

110 

244.15 

■iUI 

9x9  18 

108 

164.46 

1.27 

162.67 

9x9  19 

128 

173.27 

1.07 

171.72 

9x9  20 

130 

259.19 

mam 

257.70 

261.39 

9x9  21 

132 

mem 

199.71 

204.81 

9x9  22 

116 

165.38 

0.57 

164.62 

167.35 

9x9  23 

136 

198.46 

196.65 

202.37 

9x9  24 

120 

259.19 

258.58 

261.34 

9x9  25 

120 

178.98 

■KM 

177.40 

183.23 

9x9  26 

124 

182.78 

1.61 

180.81 

186.85 

9x9  27 

132 

199.87 

0.99 

198.22 

204.15 

9x9  28 

102 

163.11 

0.64 

162.26 

164.83 

9x9  29 

114 

216.46 

0.84 

215.14 

218.58 

9x9  30 

118 

211.37 

0.62 

210.42 

213.34 

To  determine  if  the  solution  found  using  the  simulated  annealing  algorithm  was  indeed  optimal,  the 
following  integer  programming  (IP)  model  of  SRCFP  is  formulated,  based  on  the  Miller-Tucker-Zemlin 
formulation  (Miller  et  al.  1960)  of  the  traveling  salesman  problem. 
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i&V 


Min  S  DyXy 

I  =1 

ysi 

E  =  1  ^  S 

JeL* 

Zxy=l  jeL^ 

ieL 

ZXij=l  j^L- 

iet 

p^+\-k{^-Xy^<  Pj  for  all  (/,  j)  e  A,  with  j 

Xy  +  Xj^  <  1  for  all  (i,j)  e  A,  with  i  <  j 

x^  e{0,l}  forall  (/,7)e^ 

Pie{l,...,k}  i  =  l,...,k 

where  L  =  {//,  4}  is  the  set  of  single  unit  cut  and  fill  locations,  V'  is  the  set  of  unit  cut  locations,  L' 

is  the  set  of  unit  fill  locations,  k  is  the  number  of  unit  cut  and  fill  locations,  A  is  the  set  of  arcs  from  unit 
cut  locations  to  unit  fill  locations  and  from  unit  fill  locations  to  unit  cut  locations,  pi  is  the  position  of  U  in 
the  Hamiltonian  circuit,  and  Xy  is  one  if  //  immediately  precedes  Ij  in  the  Hamiltonian  circuit,  otherwise  Xjj 
is  zero. 

ILOG  OPL  Studio  3.0  was  used  to  generate  the  integer  programs  for  the  test  problems  and  CPLEX 
6.6  was  used  to  attempt  to  solve  the  resulting  IP’s.  The  CPLEX  parameters  were  set  at  their  default 
values,  with  the  time  limit  set  at  five  hours  and  the  upper  cutoff  parameter  set  equal  to  the  best  value 
found  by  the  simulated  annealing  algorithm  plus  0.001.  The  problems  were  run  on  a  Pentium  m  550 
MHZ  processor  with  128Mb  of  random  access  memory. 

CPLEX  can  be  used  to  solve  the  randomly  generated  5x5  problem  instances.  However,  for  a 
randomly  generated  7x7  problem  instance,  CPLEX  was  terminated  after  four  hours  of  computation 
because  it  ran  out  of  random  access  memory  to  store  the  subproblems.  During  the  search,  CPLEX 
explored  nearly  34,000  subproblems,  but  did  not  find  a  better  solution  than  the  best  one  found  by  the 
simulated  annealing  algorithm.  Furthermore,  it  appears  that  the  problem  cannot  be  solved  in  any 
reasonable  length  of  time  even  on  a  faster  machine  with  more  memory,  since  there  were  approximately 
31,000  subproblems  waiting  to  be  explored  at  the  time  of  termination,  and  this  number  was  still 
increasing. 

Since  CPLEX  was  unable  to  solve  the  larger  problems,  and  thus  determine  if  the  best  solutions  found 
by  the  simulated  annealing  algorithms  are  optimal,  the  Held-Karp  1-tree  lower  bound  (see  Held  and  Karp 
1970,  1971)  was  computed  for  each  of  these  problems  in  order  to  compute  an  upper  bound  on  the 
distance  between  the  best  value  found  by  the  simulated  annealing  algorithm  and  the  optimal  value.  A 
brief  description  of  this  lower  bound  is  included  here.  Every  solution  of  SRCFP  is  a  Hamiltonian  circuit, 
which  can  be  viewed  as  a  spanning  tree  of  nodes  2,  3,,..,  k  together  with  two  edges  incident  to  node  1. 
Therefore,  the  cost  of  a  minimum  spanning  tree  of  nodes  2,  3,...,  k  plus  the  cost  of  the  two  cheapest 
edges  incident  to  node  1  is  a  lower  bound  for  the  problem.  This  lower  bound  can  be  improved  by  using 
Lagrange  multipliers  on  the  nodes  of  the  graph.  Held  et  al.  (1974)  show  that  subgradient  optimization 
can  be  used  to  converge  to  the  optimal  set  of  Lagrange  multipliers.  Tables  4.4  through  4.6  contain  the 
best  value  found  by  the  simulated  annealing  algorithm  and  the  Held-Karp  1-tree  lower  bounds  for  the 
ninety  problem  instances.  They  also  contain  the  gap  (defined  as  the  difference  between  the  best  value 
found  by  the  simulated  annealing  algorithm  and  the  lower  bound  divided  by  the  lower  bound). 
Moreover,  a  greedy  algorithm  is  presented  as  an  upper  bound  to  assess  the  effectiveness  of  simulated 
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annealing.  The  greedy  algorithm  implemented  takes  an  initial  unit  cut  location,  Cj,  and  moves  to  the 
nearest  (i.e.,  closest  in  distance)  unvisited  unit  fill  location  =  argmin{Z)(ci,  ij  =  1,  2,...,  m,fj 

has  not  been  visited}.  From  this  unit  fill  location,  the  greedy  algorithm  moves  to  the  nearest  unvisited 
cut  location  (i.e.,  c,  =  argmin{T)07,  c,),  i,j  =  1,  2,...,  m,  Cf  has  not  been  visited).  This  process  continues 
until  the  construction  project  site  is  leveled,  resulting  in  a  paired  route  of  alternating  unit  cut  and  unit  fill 
locations.  This  straightforward  technique  provides  a  route  that  although  rarely  optimal  has  an  objective 
function  value  that  is  typically  much  below  that  of  a  randomly  generated  solution. 


Table 

Held-Karp  1-Tree  Lower  Bounds 
(5  Unit  Cut  Locations  and 

4.4 

and  Greedy  Algorithm  Results 

5  Unit  Fill  Locations) 

Problem 

SA  Best  Value 

Lower  Bound 
(LB) 

Gap  between  SA 
Best  Value  and  LB 

Greedy 

Algorithm 

5x5  1 

41.14 

41.136 

0.0000 

46.54 

5x5  2 

47.37 

47.37 

0.0000 

54.19 

5x5  3 

60.97 

60.96 

0.0002 

69.83 

5x5  4 

49.00 

48.74 

0.0053 

51.36 

5x5  5 

58.80 

58.80 

0.0000 

75.78 

5x5  6 

68.02 

67.73 

0.0043 

78.48 

5x5  7 

44.50 

44.41 

0.0021 

51.94 

5x5  8 

64.47 

64.47 

0.0000 

72.37 

5x5  9 

77.24 

77.16 

0.0009 

90.94 

5x5  10 

48.55 

48.25 

0.0063 

57.14 

5x5  11 

54.65 

54.54 

0.0021 

73.60 

5x5  12 

43.66 

43.67 

0.0000 

52.65 

5x5  13 

51.02 

50.81 

0.0042 

59.84 

5x5  14 

60.78 

60.78 

0.0000 

71.75 

5x5  15 

62.39 

62.34 

0.0008 

72.39 

5x5  16 

53.90 

53.53 

0.0068 

66.72 

5x5  17 

54.14 

54.14 

0.0000 

65.06 

5x5  18 

46.34 

46.34 

0.0000 

58.38 

5x5  19 

38.33 

38.07 

0.0067 

41.75 

5x5  20 

40.60 

40.58 

0.0006 

44.95 

5x5  21 

36.96 

36.96 

0.0000 

40.49 

5x5  22 

49.98 

49.97 

0.0003 

56.17 

5x5  23 

58.19 

58.06 

0.0023 

63.22 

5x5  24 

58.09 

58.09 

0.0000 

63.63 

5x5  25 

66.84 

66.83 

0.0001 

80.03 

5x5_26 

88.01 

87.95 

0.0007 

95.71 

5x5_27 

65.36 

65.30 

0.0010 

74.98 

5x5  28 

41.38 

41.34 

0.0008 

48.43 

5x5  29 

57.16 

57.10 

0.0011 

65.60 

5x5  30 

54.61 

54.28 

0.0062 

61.51 
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Table  4.5 

-Karp  1-Tree  Lower  Bounds  and  Greedy  Algorithm  Results 
(7  Unit  Cut  Locations  and  7  Unit  Fill  Locations) 

SA  Best  Value  Lower  Bound  Gap  between  SA 


Best  Value  and  LB 


92 _ 

30 


34 _ 

85 


34 _ 

57 


2 _ 

1 


Aleorithm 


27.10 

31.37 


136.15 

148.99 


09.57 

28.58 


76.88 

83.90 


Table  4.6 

Held-Karp  1-Tree  Lower  Bounds  and  Greedy  Algorithm  Results 
_  (9  Unit  Cut  Locations  and  9  Unit  Fill  Locations) _ 


Problem 

SA  Best  Value 

Lower  Bound 
(LB) 

Gap  between  SA 
Best  Value  and  LB 

Greedy 

Algorithm 

9x9_l 

244.62 

243.86 

0.0032 

311.84 

9x9  2 

170.17 

169.26 

0.0054 

204.40 

9x9  3 

202.25 

201.04 

0.0060 

237.06 

9x9_4 

177.00 

176.31 

0.0040 

222.43 

9x9  5 

214.34 

213.27 

0.0050 

264.99 

9x9  6 

210.53 

207.80 

0.0131 

250.71 

9x9  7 

202.53 

200.21 

0.0116 

245.63 

9x9  8 

186.83 

184.31 

0.0137 

242.87 

9x9_9 

191.87 

191.62 

0.0013 

216.29 

9x9  10 

223.66 

222.67 

0.0044 

267.37 

9x9  11 

217.93 

216.29 

0.0076 

286.19 

9x9_ 12 

242.49 

241.16 

0.0055 

302.43 

9x9  13 

170.40 

168.45 

0.0115 

217.73 

9x9  14 

168.10 

167.15 

0.0057 

198.37 

9x9  15 

216.33 

214.90 

0.0067 

284.97 

9x9  16 

153.78 

151.91 

0.0123 

204.55 

9x9  17 

242.87 

241.80 

0.0044 

289.50 

9x9  18 

162.67 

161.78 

0.0055 

196.99 

9x9  19 

171.72 

170.70 

0.0060 

194.25 

9x9  20 

257.70 

255.67 

0.0079 

314.66 

9x9  21 

199.71 

198.26 

0.0073 

240.30 

9x9  22 

164.62 

161.56 

0.0189 

203.28 

9x9_23 

196.65 

195.04 

0.0083 

259.14 

9x9  24 

258.58 

258.29 

0.0011 

305.55 

9x9  25 

177.40 

175.48 

0.0109 

228.07 

9x9  26 

180.81 

179.11 

0.0095 

233.47 

9x9  27 

198.22 

196.28 

0.0099 

256.40 

9x9_28 

162.26 

162.04 

0.0013 

197.52 

9x9_29 

215.14 

214.91 

0.0010 

262.76 

9x9  30 

210.42 

210.20 

0.0011 

252.72 

The  best  solutions  found  by  the  simulated  annealing  algorithm  had  an  average  gap  of .  1 8%,  .5 1  %  and 
.7%  for  the  randomly  generated  instances  of  problems  corresponding  to  grid  sizes  5  x  5,  7  x  7,  and  9x9, 
respectively.  Moreover,  all  optimal  solutions  found  using  simulated  annealing  were  within  .68%,  1.72%, 
and  1.89%  of  the  lower  bound  for  the  randomly  generated  instances  of  problems  corresponding  to  grid 
sizes  5  X  5,  7  X  7,  and  9x9,  respectively.  Note  that  the  gap  is  measured  relative  to  a  lower  bound,  not 
the  true  optimal  solution.  It  is  possible  that  a  significant  portion  of  the  gap  is  due  to  the  lower  bound. 
The  best  solutions  obtained  by  the  simulated  annealing  algorithm  are  uniformly  better  than  those 
obtained  by  the  greedy  algorithm. 

In  the  construction  industry,  using  earthmoving  vehicles  to  transform  an  existing  terrain  into  a  final 
grade  (often  a  level  project  site)  is  a  costly  operation  and  may  represent  a  significant  portion  of  the  over¬ 
all  construction  project  budget.  Reducing  the  total  distance  traveled  by  earthmoving  vehicles  results  in 
cost  savings  in  terms  of  fuel  consumption,  time,  and  equipment  maintenance.  This  report  summarizes  the 
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SRCFP,  which  seeks  to  identify  a  route  that  minimizes  the  total  haul  distance  traveled  by  a  vehicle 
between  cut  and  fill  locations  and  returns  the  vehicle  to  its  starting  location  once  the  final  grade  is 
achieved.  Finding  a  route  that  minimizes  the  total  haul  distance  is  computationally  hard.  Simplifying  the 
SRCFP  to  a  transportation  problem  results  in  the  optimal  assignment  of  resources,  but  does  not  provide 
an  optimal  route.  Formulating  the  SRCFP  as  a  discrete  optimization  problem  and  using  local  search 
strategies  provides  near-optimal  vehicle  routes  in  a  reasonable  amount  of  time. 

There  are  several  opportunities  to  extend  the  results  described  in  this  report.  Construction  sites  and 
quarries  are  rarely  oriented  to  allow  complete  freedom  of  maneuvering  for  all  vehicle  types.  Routes  are 
limited  to  haul  roads  or  vehicle  characteristics  that  make  some  terrain  impassible  or  prohibitively 
expensive  to  traverse,  in  terms  of  fuel  and/or  equipment  wear.  One  assumption  in  the  SRCFP 
formulation  is  that  the  paths  between  cut  and  fill  location  are  direct  lines  (i.e.,  the  equipment  is  capable 
of  negotiating  any  terrain  without  penalty).  A  more  realistic  approach  is  to  incorporate  haul  roads  and 
terrain  features  into  the  distance  matrix  such  that  a  more  realistic  construction  site  can  be  modeled. 
Work  is  in  progress  to  collect  terrain  data  fi-om  actual  construction  sites  that  can  be  incorporated  into 
future  models. 

A  natural  extension  of  the  current  problem  formulation  is  to  add  penalties  to  routes  that  force  the 
equipment  to  negotiate  undesirable  terrain.  Given  the  choice  between  avoiding  a  hill  versus  going  over  a 
hill,  the  equipment  operator’s  decision  is  based  on  cost  in  terms  of  fuel,  distance  or  feasibility. 
Incorporating  penalty  functions  for  undesirable  terrain  may  provide  more  realistic  solutions  that  are  of 
greater  interest  to  the  construction  industry.  Vehicle  performance  characteristics  and  fuel  consumption 
rates  are  available  and  can  also  be  incorporated  into  the  model.  All  the  modifications  described  here  are 
meant  to  add  realism  to  the  problem.  However,  they  require  extensive  data  collection  that  may  be 
difficult  to  obtain  from  actual  construction  projects  to  support  the  resulting  model. 

The  same  technique  used  for  optimizing  cut  and  fill  patterns  can  be  applied  to  rapid  runway  repair 
for  both  military  and  nonmilitary  airfields.  Rapid  runway  repair  is  a  concern  for  civilian  airports  that 
require  maintenance  during  limited  hours  of  non-use,  and  with  military  airfields  that  require  repair  from 
incoming  artillery,  bombs  or  missiles.  In  both  cases,  time  of  repair  and  availability  of  equipment  is  the 
major  concern,  since  busy  airports  cannot  afford  to  have  a  runway  down  for  a  significant  length  of  time 
and  military  runways  must  be  fully  mission  capable  during  times  of  conflict.  The  SRCFP  can  be 
modified  to  accommodate  multiple  types  of  equipment  such  as  earthmoving  vehicles  as  well  as  concrete 
and  asphalt  transporters.  Minimizing  travel  time  by  developing  routes  for  each  type  of  vehicle  between 
material  stockpiles  to  repair  sites  is  critical  to  the  success  of  these  operations. 

Optimal  strategies  for  conducting  searches  with  limited  assets  are  also  a  natural  extension  of  this 
model.  Different  platforms  (e.g.,  helicopters  versus  fixed  wing  or  satellite)  with  varying  degrees  of 
effectiveness  are  used  for  search  and  rescue  operations.  Due  to  the  timeliness  required  in  search  and 
rescue  operations,  efficient  use  of  available  assets  is  critical  to  success.  The  optimal  routes  developed  by 
approaching  the  SRCFP  with  local  search  strategies  can  be  extended  to  optimize  a  fleet  of  search 
platforms.  The  same  concepts  can  be  extended  to  optimizing  surveillance  assets  to  maximize  coverage  of 
a  given  location. 

5.  Other  Research  Results 

In  addition  to  the  results  reported  above,,  several  other  results  were  obtained  during  the  period  of  the 
grant.  These  results  are  briefly  discussed  here. 

Jacobson,  Kobza  and  Easterling  (2001)  present  a  new  approach  to  addressing  the  difficult  tradeoff 
between  false  alarms  and  false  clears  in  aviation  security  systems.  By  modeling  the  problem  as  a  discrete 
optimization  problem,  the  paper  establishes  that  the  resulting  problem  to  be  NP-hard,  and  then  develops 
heuristic  procedures  to  address  the  problem.  The  results  in  this  paper  have  the  potential  to  change  the 
way  in  which  aviation  security  system  information  can  be  synthesized  and  interpreted.  These  results  may 
also  have  important  applications  in  health  care  (in  diagnosing  certain  types  of  cancer)  and  information 
systems  (in  detecting  when  a  system  has  crashed).  Note  that  this  paper  also  received  the  2003  Best  Paper 
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award  in  HE  Transactions  focused  issue  on  Operations  Engineering.  Jacobson,  Bowman,  and  Kobza 
(2001)  introduce  and  analyzes  the  flight  segment  baggage  value  (FSBV)  and  the  passenger  segment 
baggage  value  (PSBV)  performance  measures  for  baggage  screening  security  systems.  The  paper  also 
includes  a  real-world  example  using  actual  flight  data  from  the  Official  Airline  Guide  (OAG)  to  illustrate 
an  application  of  the  models  and  results  presented  in  the  paper.  Virta  et  al.  (2002)  look  at  an  important 
extension  of  how  the  origination  of  selectees  can  impact  the  overall  security  of  an  air  system.  These 
results  are  particularly  noteworthy  since  two  of  the  19  hijackers  during  the  9/11/2001  terrorist  activities 
boarded  planes  in  Portland,  Maine  and  connected  to  their  flights  in  Boston.  The  results  in  this  paper 
identify  how  such  strategies  are  a  major  weakness  of  the  air  security  system,  and  quantify  the  risks 
associated  with  such  security  holes.  Jacobson  et  al.  (2003)  describe  how  discrete  optimization  models 
can  be  used  to  address  aviation  security  system  deployment  and  utilization  questions,  based  on  three 
performance  measures  that  quantify  the  effectiveness  of  airport  baggage  screening  security  device 
systems.  These  models  are  used  to  solve  for  optimal  airport  baggage  screening  security  device 
deployments  considering  the  number  of  passengers  on  a  set  of  flights  who  have  not  been  cleared  using  a 
security  risk  assessment  system  in  use  by  the  Federal  Aviation  Administration  (i.e.,  passengers  whose 
baggage  is  subjected  to  screening),  the  number  of  flights  in  this  set,  and  the  size  of  the  aircraft  for  such 
flights.  Several  examples  are  provided  to  illustrate  these  results,  including  an  example  that  uses  data 
available  from  the  Offieial  Airline  Guide.  Virta  et  al.  (2003)  consider  the  cost  and  benefits  of  various 
checked  baggage  screening  strategies.  Determining  how  to  effectively  operate  security  devices  is  as 
important  to  overall  system  performance  as  developing  more  sensitive  security  devices.  In  light  of  recent 
federal  mandates  for  100%  screening  of  all  checked  baggage,  this  paper  studies  the  tradeoffs  between 
screening  only  selectee  cheeked  baggage  and  screening  both  seleetee  and  non-selectee  checked  baggage 
for  a  single  baggage  screening  security  device  deployed  at  an  airport.  This  tradeoff  is  represented  using  a 
cost  model  that  incorporates  the  eost  of  the  baggage  screening  security  device,  the  volume  of  checked 
baggage  processed  through  the  device,  and  the  outcomes  that  occur  when  the  device  is  used.  The  cost 
model  captures  the  cost  of  deploying,  maintaining,  and  operating  a  single  baggage  screening  security 
device  over  a  one-year  period.  The  study  concludes  that  as  excess  baggage  screening  capacity  is  used  to 
screen  non-selectee  checked  bags,  the  expected  annual  cost  increases,  the  expected  annual  cost  per 
checked  bag  screened  decreases,  and  the  expected  annual  cost  per  expected  number  of  threats  detected  in 
the  checked  bags  screened  increases.  These  results  indicate  that  the  marginal  increase  in  security  per 
dollar  spent  is  significantly  lower  when  non-selectee  checked  bags  are  screened  than  when  only  selectee 
checked  bags  are  screened. 

Vazquez-Abad  and  Jacobson  (2001)  present  a  new  gradient  estimator  for  the  steady-state  expected 
sojourn  (system)  time  in  a  nonpreemptive  priority  queueing  system.  The  estimator  uses  the  concept  of  a 
phantom  system,  together  with  the  basic  ideas  in  harmonic  gradient  estimation,  to  develop  a  single 
simulation  run  estimator,  termed  the  phantom  harmonic  gradient  estimator.  The  estimator  is  shown  to  be 
strongly  consistent  and  strongly  eonsistent  in  the  average  sense  as  the  sample  size  grows.  An  upper 
bound  for  the  variance  of  the  PHG  estimator  is  presented.  This  bound  is  used  to  show  that  under  mild 
conditions,  the  variance  of  the  PHG  estimator  tends  to  zero  as  both  the  number  of  phantom  systems  and 
the  sample  size  approach  infinity.  A  variance  reduction  technique  that  simultaneously  uses  both  common 
and  antithetic  random  numbers  is  presented.  Computational  results  on  several  non-preemptive  queueing 
systems  illustrate  the  effectiveness  of  the  method. 

Sullivan  and  Jacobson  (2001)  present  new  necessary  and  sufficient  convergence  conditions  for 
generalized  hill  climbing  algorithms.  These  conditions  are  then  contrasted  with  similar  conditions  for 
simulated  annealing  algorithms,  a  particular  type  of  generalized  hill  climbing  algorithms,  to  show  that 
they  reduce  to  the  most  often  cited  and  used  conditions  for  such  algorithms.  These  results  provide 
important  insights  into  how  the  restrictive  algorithmic  constraints  imposed  by  simulated  annealing  can  be 
circumvented  using  relative  rates  at  which  global  and  local  optima  can  be  accessed  in  the  solution  space 
by  generalized  hill  climbing  algorithms. 
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Jacobson  and  Yucesan  (2002)  identify  and  investigate  common  links  between  discrete-event 
simulation  and  discrete  optimization  algorithms.  The  primary  contribution  of  the  paper  is  two  search 
problem  formulations  that  are  proven  to  be  NP-hard,  and  which  have  implication  both  in  the  design  and 
analysis  of  discrete  event  simulation  models  and  in  the  design  and  implementation  of  discrete 
optimization  algorithms.  Fleischer  and  Jacobson  (2002)  show  how  the  Boltzmann  distribution  used  in 
the  steady-state  analysis  of  the  simulated  annealing  algorithm  gives  rise  to  several  scale  invariant 
properties.  Scale  invariance  is  first  presented  in  the  context  of  parallel  independent  processors  and  then 
extended  to  an  abstraet  form  based  on  lumping  states  together  to  form  new  aggregate  states.  These 
lumped  or  aggregate  states  possess  all  of  the  mathematical  characteristics,  forms  and  relationships  of 
states  (solutions)  in  the  original  problem  in  both  first  and  second  moments.  These  scale  invariance 
properties  therefore  permit  new  ways  of  relating  objective  function  values,  conditional  expectation 
values,  stationary  probabilities,  rates  of  change  of  stationary  probabilities  and  eonditional  variances. 
Such  properties  therefore  provide  potential  applications  in  analysis,  statistical  inference  and 
optimization. 

Sewell  and  Jacobson  (2003)  is  part  of  an  on  going  research  activity  into  the  application  of  operation 
research  tools  in  the  health  eare  domain.  The  Recommended  Childhood  Immunization  Schedule  has 
become  sufficiently  crowded  that  the  prospect  of  adding  additional  vaccines  to  this  schedule  may  not  be 
well  received  by  either  health-  care  providers  or  parents/guardians.  This  has  encouraged  vaccine 
manufacturers  to  develop  combination  vaccines  that  can  permit  new  vaccines  to  be  added  to  the  schedule 
without  requiring  children  to  be  exposed  to  an  unacceptable  number  of  injections  during  a  single  clinic 
visit.  This  paper  develops  an  integer  programming  model  to  assess  the  economic  premium  that  exists  in 
having  combination  vaccines  available.  The  results  of  this  study  suggest  that  combination  vaccines 
provide  a  cost  effective  alternative  to  individual  vaccines  and  that  further  developments  and  innovations 
in  this  area  by  vaceine  manufacturers  can  provide  significant  economic  and  societal  benefits.  Jacobson, 
Kamani,  and  Sewell  (2003)  report  the  results  of  reverse  engineering  a  vaccine  selection  algorithm  to 
evaluate  the  economic  value  of  a  hepatitis  B  -  Haemophilus  influenzae  type  B  combination  vaccine  that 
is  currently  under  federal  contract  in  the  United  States.  This  analysis  captures  the  tradeoff  between  the 
cost  assigned  to  administering  an  injection  and  the  price  of  the  vaccine  that  earns  it  a  place  in  the  lowest 
overall  cost  formulary.  Jacobson  and  Sewell  (2002)  report  other  such  results  using  Monte  Carlo 
simulation. 

Swisher  et  al.  (2003)  present  a  survey  of  the  literature  for  two  widely  used  classes  of  statistical 
methods  for  selecting  the  best  design  fi'om  among  a  finite  set  of  k  alternatives:  ranking  and  selection 
(R&S)  and  multiple  comparison  procedures  (MCPs).  A  comprehensive  survey  of  each  topic  is  presented 
along  with  a  summary  of  recent  unified  R&S-MCP  approaches.  Procedures  are  recommended  based  on 
their  statistical  efficiency  and  ease  of  application;  guidelines  for  procedure  application  are  offered. 

Armstrong  and  Jacobson  (2003)  examine  the  complexity  of  global  verification  for  MAX-SAT, 
MAX-A:-SAT  (for  k  >  3),  Vertex  Cover,  and  the  Traveling  Salesman  Problem.  These  results  are  obtained 
by  adaptations  of  the  transformations  that  prove  such  problems  to  be  NP-complete.  The  class  of 
problems  PGS  is  defined  to  be  those  discrete  optimization  problems  for  which  there  exists  a  polynomial 
time  algorithm  sueh  that  given  any  solution  co,  either  a  solution  can  be  found  with  a  better  objective 
function  value  or  it  can  be  concluded  that  no  such  solution  exists  and  co  is  a  global  optimum.  This  paper 
demonstrates  that  if  any  one  of  MAX-SAT,  MAX-A:-SAT  (for  k  ^  3),  Vertex  Cover,  or  Traveling 
Salesman  Problem  are  in  PGS,  then  P  =  NP. 
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